Vincent Lemort
email: vincent.lemort@ulg.ac.be

Sylvain Quoilin
email: squoilin@ulg.ac.be

Cristian Cuevas
email: crcuevas@udec.cl

Jean Lebrun
email: j.lebrun@ulg.ac.be

Testing and modeling a scroll expander integrated into an organic rankine cycle

Keywords: Scroll, expander, Rankine cycle, heat recovery

Organic Rankine Cycles (ORC's) are particularly suitable for recovering energy from low-grade heat sources. This paper first presents the results of an experimental study carried out on a prototype of an open-drive oil-free scroll expander integrated into an ORC working with refrigerant HCFC-123. By exploiting the overall expander performance measurements, the eight parameters of a scroll expander semi-empirical model are then identified. The model is able to compute variables of first importance such as the mass flow rate, the delivered shaft power and the discharge temperature, and secondary variables such as the supply heating-up, the exhaust cooling-down, the ambient losses, the internal leakage and the mechanical losses. The maximum deviation between the predictions by the model and the measurements is 2% for the mass flow rate, 5% for the shaft power and 3K for the discharge temperature. The validated model of the expander is finally used to quantify the different losses and to indicate how the design of the expander might be altered to achieve better performances. This analysis pointed out that the internal leakages and, to a lesser extent, the supply pressure drop and the mechanical losses are the main losses affecting the performance of the expander.

ACCEPTED MANUSCRIPT

Among the solutions currently proposed, the Organic Rankine Cycle (ORC) is the most widely used. Its main advantages are its simplicity and its commonly available components. In such a cycle, the working fluid is an organic product that presents a lower ebullition temperature than water, allowing reduced evaporating temperatures.

Unlike traditional power cycles, local and small-scale energy production is made possible with this technology (Kane [START_REF] Kane | Small hybrid solar power system[END_REF], Quoilin [START_REF] Quoilin | Experimental study and modeling of a low temperature Rankine Cycle for small scale cogeneration[END_REF]).

The expander is a key element of the ORC. The choice of the expander strongly depends on the operating conditions and on the size of the facility. Two main types of machines can be distinguished: the dynamic (turbo) and displacement (volumetric) type.

Displacement type machines are more appropriate to the small-scale ORC unit under investigation in this paper, because they are characterized by lower flow rates, higher pressure ratios and much lower rotational speeds than turbo-machines (Persson [3]).

Moreover, these machines can tolerate two-phase conditions, which may appear at the end of the expansion in some operating conditions.

Among positive displacement machines, the scroll machine is a good candidate for the ORC application, because of its reduced number of moving parts, reliability, wide output power range, and broad availability (Zanelli and Favrat [START_REF] Zanelli | Experimental Investigation of a Hermetic Scroll Expander-Generator[END_REF]). Moreover, it is a proven technology in compressor mode due to its extensive use in refrigeration and airconditioning industry.

However, up to now, the use of scroll machines in expander mode has mainly been limited to experimental work and so far numerous scroll expander prototypes have been tested for different fluids.

For instance, Yanagisawa et al. [START_REF] Yanagisawa | Performance of an oil-free scrolltype air expander[END_REF] carried out an experimental study on an oil-free scroll-type air expander. They observed that the performance is heavily reduced by the mechanical loss, but leakage loss becomes significant as the rotational speed decreases.

Maximal achieved volumetric and isentropic effectiveness was respectively 76% and 60%. Zanelli and Favrat [START_REF] Zanelli | Experimental Investigation of a Hermetic Scroll Expander-Generator[END_REF] carried out an experimental investigation on a hermetic scroll expander-generator fed with refrigerant R134a. The machine produced a power ranging from 1.0 to 3.5 kW with a maximal isentropic effectiveness of 65%. Kane [START_REF] Kane | Small hybrid solar power system[END_REF] designed, built and tested a prototype of hybrid solar thermal power plant associating solar collectors, cogeneration engines and two superposed Organic Rankine Cycles, equipped with hermetic scroll expanders. Manzagol et al. [START_REF] Manzagol | Cryogenic scroll expander for Claude cycle with cooling power of 10 to 100 watts at 4.2 K[END_REF] studied a cryogenic scroll expander used for a 10 L/h helium liquefier. The expander was tested on a Brayton cycle refrigerator and reached an isentropic effectiveness of 50 to 60% for supply gas conditions of 35 K and 0.7 MPa. Xiaojun et al. [START_REF] Xiaojun | Research on a Scroll Expander Used for Recovering Work in a Fuel Cell[END_REF] investigated the possibility to recover work in a fuel cell by means of a scroll expander. The work produced by the expander could be provided to the compressor. Hugenroth et al. [START_REF] Hugenroth | Experimental investigation of a liquidflooded Ericsson cycle cooler[END_REF] developed a novel approach to implementing a gas Ericsson cycle heat pump, using liquid flooding of the compressor and of the expander to approach isothermal compression and expansion process. Scroll machines were selected due to their ability to compress and expand a mixture of oil and gas. Aoun and Clodic [START_REF] Aoun | Theoretical and experimental study of an oil-free scroll type vapor expander[END_REF] investigated the same expander as Yanagisawa et al. [START_REF] Yanagisawa | Performance of an oil-free scrolltype air expander[END_REF], but operating with steam. In order to improve the volumetric performance, original tip seals were replaced by PTFE tip seals, better adapted for high temperature applications and presenting lubricating properties. Maximal achieved volumetric and isentropic effectiveness was respectively 63% and 48%.

This paper presents a semi-empirical simulation model of a scroll expander. Unlike deterministic models (that require the exact knowledge of the geometry of the machine), the model proposed only involves a limited number of parameters. By using experimental data these parameters will be identified in the case of a prototype of an expander fed with HCFC-123.

EXPERIMENTAL INVESTIGATION

Description of the test rig

An experimental investigation is carried out on a scroll expander integrated into a prototype of an ORC system working with refrigerant HCFC-123. This paper only focuses on the experimental characterization of the scroll expander and not of the entire ORC system.

The test rig is schematically described in the expander performance is given in Table 1.

Description of the tests

The experimental data covers a large range of operating conditions. In total, 39 points are achieved. Table 2 indicates what are the imposed and the obtained variables during 2 indicates that the vapor is superheated at the expander supply for all of the performance points.

Results

The ranges of the main measured variables are given in Table 2. The maximum delivered shaft power is 1.82 kW and the maximum achieved overall isentropic ACCEPTED MANUSCRIPT 10 effectiveness is 68%. This overall effectiveness is defined by the ratio of the measured shaft power and the isentropic power (Eq. ( 1)). The latter is the product of the measured flow rate by the expansion work associated to an isentropic expansion from supply conditions to the exhaust pressure:

) h h .( M W w . M W W W s ,
- = = = ε (1)
The evolution of the overall isentropic effectiveness with the pressure ratio imposed to the expander for the three rotational speeds is given in Figure 2. The sharp decrease in effectiveness for the low pressure ratios is due to over-expansion losses. Error bars are associated with calculated uncertainties. Provided measurements are uncorrelated and random, the uncertainty U Y on the variable Y is calculated as function of the uncertainties U Xi on each measured variables X i by Eq. (2) (Klein [START_REF] Klein | Engineering Equation Solver, F-Chart Software[END_REF]).

∂ ∂ = i X i Y i U X Y U 2 2 (2)
Figure 2 also shows that the maximum achieved pressure ratio is around 5.5. Further experimental work should allow evaluating the performance of the expander for higher pressure ratios.

The volumetric performance of the expander is represented by the filling factor (Zanelli and Favrat [START_REF] Zanelli | Experimental Investigation of a Hermetic Scroll Expander-Generator[END_REF]). The latter is defined as the ratio between the measured mass flow rate and the mass flow rate theoretically displaced by the expander (Eq. ( 3)). The filling factor increases with internal leakage and supply cooling down, but decreases with the supply pressure drop.

s su meas meas V v . M = φ (3)
The evolution of the filling factor with the expander supply pressure for the three rotational speeds is given in Figure 3. The filling factor could be expected to increase with the supply pressure due to larger internal leakage flows. However, this trend is balanced by the simultaneous effect of the supply pressure drop. It can be observed that the lower the rotational speed, the higher the filling factor, due to the larger relative impact of the internal leakage and due to smaller supply pressure drop.

EXPANDER MODEL

The semi-empirical model of a scroll expander described hereunder is deduced from the one proposed by Winandy et al. [START_REF] Winandy | Experimental analysis and simplified modelling of a hermetic scroll refrigeration compressor[END_REF] for hermetic scroll compressors. It is similar to the model of Kane [START_REF] Kane | Small hybrid solar power system[END_REF], but describes internal and external heat transfers (and requires three From the description given in Figure 4, it can be observed that the heat transfers, the supply pressure drop and the internal leakage are fictitiously dissociated from the actual expansion process (su,2 ex,2). All the processes mentioned hereunder are described in the following paragraphs.

Supply pressure drop

The supply pressure drop ( ) accounts for all pressure losses encountered by the fluid from the suction line to the suction chamber. Major pressure losses are associated to two following phenomena (Yanagisawa et al. [START_REF] Yanagisawa | Performance of an oil-free scrolltype air expander[END_REF]):

-As shown in Figure 5, during part of the suction process, the expander suction port is blocked by the tip of the orbiting scroll, reducing the effective suction port.

-At the end of the suction process, the flow passage between the central portion of the suction chamber and the two adjacent crescent-shaped portions is progressively reduced to zero. However, the tip seal does not extend to the end of the scrolls, which increases the flow between the central and the adjacent chambers and attenuates the supply pressure drop.

The lumped supply pressure drop is computed by comparison to the isentropic flow through a converging nozzle, whose cross sectional area A su is a parameter to identify.

Because of the steady-state nature of the model, this cross sectional area represents an average value of the suction port effective area over the entire suction process (that extends over one shaft revolution).

By combination of the equations of mass and energy conservation through the nozzle, the mass flow rate entering the expander can be expressed as:

( ) 1 2 , su su su , thr su h h v A M - = (4)
Both the specific enthalpy and the specific volume at the nozzle throat are function of the pressure P su,1 at the throat and of the supply specific entropy s su . Knowing the mass flow rate entering the expander, Eq. ( 4) can be used to compute the pressure P su,1 .

Supply and exhaust heat transfers

The main heat transfer mechanisms inside the scroll expander occur between: 1) the expander shell and the fluid in the supply and exhaust pipes; 2) the scrolls (fixed and orbiting) and the fluid in the suction, expansion and discharge chambers; 3) between the shell and the ambient.

Both supply and exhaust heat transfers are computed by introducing a fictitious metal envelope of uniform temperature T w . This fictitious envelope represents the metal mass associated to the expander shell, the fixed and the orbiting scrolls. The supply heat transfer is given by:

( ) ( ) w , su p c . M AU , su , su su T T c M e h h M Q p su - ⋅ ⋅ ⋅ - = - ⋅ = - 1 2 1 1 (5)
The supply heat transfer coefficient AU su is assumed to vary with the mass flow rate according to:

8 . 0 , . = n n su su M M AU AU ( 6
)
where AU su,n is the nominal heat transfer coefficient corresponding to the nominal mass flow rate n M . This relationship can be justified by the Reynold's analogy for a turbulent flow through a pipe (Incropera and DeWitt [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF]) by assuming that the fluid properties, not included in this expression, remain unchanged. Exhaust heat transfer is described in a similar way.

Internal leakage

There are two different leakage paths in a scroll compressor/expander: the radial leakage is due to a gap between the bottom or the top plate and the scrolls and the flank leakage results from a gap between the flanks of the scrolls (Halm, [START_REF] Halm | Mathematical Modeling of Scroll Compressors[END_REF]). In this modeling, all the leakage paths are lumped into one unique fictitious leakage clearance, whose cross-sectional area A leak is a parameter to identify. The leakage flow rate can be computed by reference to the isentropic flow through a simply convergent nozzle, whose throat area is A leak . The pressure at the inlet of the nozzle is P su,2 . The throat pressure corresponds to the maximum between exhaust and critical pressures:

( )

leak crit ex leak thr P P MAX P , 2 , , , = (7)
The critical pressure P crit,leak is computed by considering the refrigerant vapor as a perfect gas:

+ = -1 2 , ,
M r V . v N M v V . N M v V M M M - = - = - = - = 2 2 2
(10)

Internal expansion

One working cycle of the scroll expander includes three processes: suction, expansion and discharge (Halm,[START_REF] Halm | Mathematical Modeling of Scroll Compressors[END_REF]). During the suction process, the suction chamber is in communication with the suction line and the fluid flows into the chamber. The expansion process is initiated when the suction chamber ceases to be in communication with the suction line. The discharge process begins when the discharge chambers enter in communication with the discharge line.

The internal power in W produced by the machine is the summation of the suction ( suc W ), the expansion ( expan W ) and discharge ( dis W ) powers. Moreover, since the internal expansion process (su,2 ex,2) is assumed adiabatic, the internal power can be expressed as function of initial and final enthalpy states:

) h h ( M W W W W , ex , su in dis suc in 2 2 expan - = + + = (11)

Suction power

The energy balance across the suction chamber, between the beginning and the end of the suction process yields:

exp , s , su , su , su in , su , su in suc V P v P M ) u h ( M W 2 2 2 2 2 = = - = (12)

Expansion power

The expansion process is assumed to be adiabatic and reversible. Writing the energy balance between the beginning and the end of this process yields:

( )

ad ad ad su su su in v P h v P h M W + - - = 2 2 2 expan (13)

Under-and over-expansion

Under-expansion occurs when the internal pressure ratio imposed by the expander (P su2 /P ad ) is lower than the system pressure ratio (P su2 /P ex2 ). In that case, the pressure in the expansion chambers at the end of the expansion process (P ad ) is higher than the pressure in the discharge line. This is represented in Figure 6 (a). The modeling assumes that there is no pressure drop through the discharge port. In order to equalize the pressures in the discharge chambers (the former expansion chambers) and in the discharge line, some fluid ( in M ∆ in Eq. ( 14)) has to flow out of the discharge chambers.

The modeling assumes it is achieved instantaneously as soon as the expansion chambers open onto the discharge line. The energy balance over the discharge chamber can be expressed as follows:

( )

2 2 , ex in ad in , ex in in h M u M u M M ∆ - = - ∆ - ( 14
)
Over-expansion occurs when the internal pressure ratio imposed by the expander is higher than the system pressure ratio (Figure 6 (b)). The energy balance over the discharge chamber is then similar to the one already established for the under-expansion process:

( )

2 2 , ex in ad in , ex in in h M u M u M M ∆ = - ∆ + ( 15
)
There is no work directly associated to the under-and over-expansion throttling.

Discharge power

As for the suction process, the equation of conservation of energy is applied on the discharge chamber between the beginning and the end of the discharge process, yielding:

( ) exp s, 2 2 2 V P r v P M M W , ex in , v , ex , ex in in dis - = ∆ ± - = (16)

Internal power

The internal power can be found by combining Equations ( 12), ( 13) and ( 16):

( ) ( ) [ ] 2 2 in , ex ad ad ad , su in P P v h h M W - + - = (17)

Shaft power

Mechanical losses loss W are due to friction between the scrolls and losses in the bearings. In the present modeling, all these losses are lumped into one unique mechanical loss torque T loss , that is a parameter to identify. The modeling assumes that this torque is independent of the rotational speed (Yanagisawa et al. [START_REF] Yanagisawa | Performance of an oil-free scrolltype air expander[END_REF]). Accordingly, the expander shaft power can be defined by: ( )

loss rot in loss in loss in sh T . N . . W W W W W W π 2 - = - = - = (
amb w amb amb T T AU Q - = (19)
Mechanical losses are directly "injected" into the envelope. In steady-state regime, the envelope temperature is computed by expressing its energy balance:

0 = - + - amb su ex loss Q Q Q W (20)

VALIDATION OF THE EXPANDER MODEL

The validation of the expander model is achieved in two steps. First, the values of the parameters of the model are identified on the basis of the measurements. Then, predictions by the model are compared to measurements. Relative deviations are then used to determine the quality of fit.

Identification of the parameters of the model

The parameter identification process is illustrated in the flow chart given in Figure 7.

The input variables of the model are the supply pressure, the supply temperature, the exhaust pressure and the rotational speed of the expander. The model calculates the mass flow rate displaced by the expander, the delivered mechanical power and the exhaust temperature. The model only necessitates eight parameters, which are summarized in Table 3. The nominal mass flow rate n M is only introduced as a reference to define the nominal heat transfer coefficients AU su,n and AU ex,n . Imposing the supply pressure as an input variable and the mass flow rate as an output variable is purely a convention. In fact, the mass flow rate could be imposed as an input (in the experimental set-up, it is actually imposed by the pump) and the supply pressure would be predicted by the model.

The parameters of the model are identified by minimizing a global error function accounting for the errors on the prediction of the mass flow rate, shaft power and exhaust temperature (which are the main output variables of the model): This minimization process is carried out by means of a genetic algorithm available in the EES software (Klein [12]). The identified parameters are listed in Table 3. The evolutions of the shaft power, measured and predicted by the model, with the pressure ratio are compared in Figure 9. Here also, it can be observed that the agreement is good. The maximum deviation between the model predictions and the measurements is 5%. The model predicts the exhaust temperature within 3 K. However, the model seems to slightly overestimate the exhaust temperature for high expander mean temperatures and to slightly underestimate it for the lower mean temperatures. This figure also shows the evolution of the predicted exhaust temperature if the model did not account for ambient losses. The deviation between the predicted exhaust temperature and the measured one increases with the mean temperature. This confirms that the model should account for the ambient losses to better predict the exhaust temperature.

+ + = - - - - ( ) ( ) ( )

Validation of the model

MODEL ANALYSIS

The validated model of the expander is used to quantify the different losses and to indicate how the design of the expander might be altered to achieve better performances.

Figure 11 shows the evolution of the overall isentropic effectiveness with the pressure ratio imposed to the expander. Operating conditions, related to one of the measured performance points, correspond to a supply pressure of 10.03 bar, an exhaust pressure of 2.01 bar, a supply temperature of 142°C and a rotational speed of 2296 rpm. The measured overall isentropic effectiveness for this operating point is also indicated in Figure 11 (with error bars associated with measurement uncertainties).

The evolution at the top of the figure is predicted by a model that only accounts for under-and over-expansion losses. The effectiveness goes through 1.0 at a pressure ratio equal to the internal pressure ratio (P ad =P ex ). For smaller and larger pressure ratios, the fluid is over-expanded and under-expanded respectively. Experimental results showed that the maximum achieved pressure ratio was around 5.5. If much larger pressure ratios were imposed, a machine with a larger built-in volume ratio than 4.05 would yield better performances.

The overall isentropic effectiveness decreases when accounting for the heat transfers in the modeling because of the supply cooling down of the fluid. This decrease could slightly be tempered by better insulating the expander (in the limiting case, AU amb =0 W/K).

Introducing the mechanical losses and the supply pressure drop in the modeling largely reduces the isentropic effectiveness. Further work should investigate the possibility of reducing mechanical losses by using better adapted tip seals. As already mentioned by [START_REF] Yanagisawa | Performance of an oil-free scrolltype air expander[END_REF], the supply pressure drop is an inherent characteristic of the scroll machine. A detailed modeling of the expander should answer the question of how to reduce this pressure drop by modifying the expander geometry.

The internal leakage is responsible for the major part of the performance loss. Under the assumption that the tip seals work correctly (they seal the radial gap between the tip of each scroll and the plate of the opposite scroll), the identified leakage area (Table 3) may be explained by a large flank clearance between the two scrolls. This large clearance is characteristic of scroll machines operating with a kinematically rigid configuration, where a flank gap is maintained.

CONCLUSIONS

An experimental study was carried out on a prototype of a scroll expander integrated into an Organic Rankine Cycle. The evolution of the expander performance, expressed in terms of global isentropic effectiveness and filling factor, with the operating conditions has been investigated. The tested prototype achieved a maximum isentropic effectiveness of 68%.

Based on experimental results, the eight parameters of a semi-empirical simulation model of the expander have been identified. The model, with its parameters identified, predicts the mass flow rate, the shaft power and the exhaust temperature with a good accuracy.

The validated model of the expander was used to quantify the different losses and to indicate how the design of the expander might be altered to achieve better performances. This analysis pointed out that the internal leakages and, to a lesser extent, the supply pressure drop and the mechanical losses are the main losses affecting the performance of the expander.

Due to its accuracy, low computational time and robustness, the model can easily be integrated into a Rankine cycle simulation model. The variation of the performance of the cycle with the performance of the expander could therefore be investigated. The model could also be used as a pre-design tool for estimating the main characteristics of the expander.

  Due to environmental constraints, the interest in low grade heat recovery has grown dramatically in the past decades. An important number of new solutions have been proposed to generate electricity from low temperature heat sources and are now applied to very diverse fields such as solar thermal power, geothermal heat sources, engine exhaust gases, domestic boilers, etc.

Figure 1 .

 1 The expander is originally an opendrive oil-free air scroll compressor. The only modifications made to the original compressor were to remove the cooling fan and to insulate the machine. The scroll expander has a kinematically rigid configuration, which maintains a small clearance gap between the scroll flanks (Peterson et al. [10]). Moving tip seals are embedded in a groove at the tip of both scrolls to reduce radial leakages (Inaba et al. [11]). The expander drives an asynchronous machine through two belt-and-pulley couplings and a torque-meter. Using an asynchronous machine is a convenient way to impose the rotational speed of the expander. The boiler of the cycle is made up of three plate heat exchangers HX1, HX2 and HX3 in series, fed with hot air. The condenser is made up of two plate heat exchangers in parallel and is fed with chilled water. A diaphragm pump drives the refrigerant through the cycle. Its displacement can be adjusted, which allows controlling the refrigerant flow rate through the cycle. The expander mechanical power is determined by measuring simultaneously the rotational speed and the torque developed at the torque-meter shaft. The transmission belt efficiency is estimated. The refrigerant flow rate is measured by a Coriolis flow meter located downstream of the pump. Temperatures are measured by means of T-type thermocouples (copper-constantan) and pressures are measured with piezoresistive pressure transmitters. The accuracy of the measurements involved in the evaluation of

  the tests. The refrigerant flow rate flowing through the expander is imposed by varying the pump displacement. The expander supply temperature is controlled by modifying the air flow rates and temperatures in the boiler. The expander rotational speed is set to 3 different values (1771, 2296, 2660 rpm) by modifying the diameter of the pulley at the asynchronous machine shaft. The expander exhaust pressure is imposed by adjusting the water flow rate through the condenser. Table

  additional parameters). The model has already been partially validated by tests with steam as working fluid (Lemort et al. [14]). In this paper, the model parameters are identified for the expander under investigation, integrated into a Rankine cycle and fed with HCFC-123. The conceptual scheme of the expander model is shown in Figure 4. In this model, the evolution of the fluid through the expander is decomposed into the following consecutive steps: a) Adiabatic supply pressure drop (su su,1), b) Isobaric supply cooling-down (su,1 su,2), c) Adiabatic and reversible expansion to the "adapted" pressure imposed by the built-in volume ratio of the machine (su,2 ad), d) Adiabatic expansion at a constant machine volume (ad ex,2), e) Adiabatic mixing between supply and leakage flows (ex,2 ex,1), f) Isobaric exhaust cooling-down or heating-up (ex,1 ex).

3 . 4 .

 34 supply pressure drop, the mass and energy conservation equations through the nozzle are combined to express the leakage mass flow rate: Displaced mass flow rate As shown in Eq. (10), the internal mass flow rate in M is the difference between the mass flow rate M entering the expander and the leakage mass flow rate leak M . The entering mass flow rate is the volume flow rate exp , s V divided by the specific volume of the fluid v su,2 after pressure drop and cooling down. The volume flow rate is the swept volume exp , s V multiplied by the expander rotational speed N. The swept volume in expander mode is equal to the one in compressor mode exp , s V divided by the built-in volume ratio of the machine r v,in .

18) 3 . 7 .

 37 Heat balance over the expanderACCEPTED MANUSCRIPT 20The ambient losses are computed by introducing a global heat transfer coefficient AU amb between the envelope and the ambient:

Figure 8

 8 Figure 8 compares the evolutions of the mass flow rate measured and predicted by the

Figure 10

 10 Figure 10 compares the evolutions of the exhaust temperature (measured and predicted

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 123456 Figure 1: Schematic representation of the ORC test bench

Figure 7 :Figure 8 :

 78 Figure 7: Flow chart of the parameter identification process

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Evolution of the expander shaft power (measured and predicted by the model) with the imposed pressure ratio

Table 1 : Accuracy of the measurement devices

 1

	Variable			Device uncertainty
	Torque	T tqm	± 0.1 N.m
	Temperature	T	± 0.3 K
	Supply pressure	P su	± 0.05 bar
	Exhaust pressure	P ex	± 0.025 bar
	Mass flow rate	M	meas	± 0.1%

Table 2 : Ranges of the main measured variables

 2

	Imposed					Obtained
	M	meas	T su,meas	T su,meas P ex,meas N	W	sh	meas ,	P su,meas T ex,meas s,meas	meas
	[g/s]	[°C]	[K]	[bar]	[min -1 ] [W]	[bar]	[°C]	[-]	[-]
	45.13 -	101.7 -	2.2 -	1.38 -	1771 -	382 -	5.45 -	66.4 -	0.424 -	1.067 -
	86.25	165.2	68.9	2.66	2660	1820	11.12	128.0	0.680	1.336

Table 3 : Identified parameters of the expander model

 3

	Heat transfer coefficient with the ambient AU amb	6.4 W/K
	Supply heat transfer coefficient	AU su,n	21.2 W/K
	Exhaust heat transfer coefficient	AU ex,n	34.2 W/K
	Nominal mass flow rate	M	n	0.12 kg/s
	Leakage area	A leak	4.6 mm 2


Mohamed Rady
email: makrady@yahoo.com

GRANULAR PHASE CHANGE MATERIALS FOR THERMAL ENERGY STORAGE: EXPERIMENTS AND NUMERICAL SIMULATIONS

The present article reports on the utilization of granular phase change composites (GPCC) of small particle diameter (1-3 mm) in latent heat thermal energy storage (LHTES) systems. The phase changing parameters (phase change temperature, latent heat, and energy storage capacity) of GPCC have been determined using differential scanning calorimeter (DSC) and temperature-history methods. Further analysis of measurement results has been conducted to describe the evolution of latent heat with temperature during phase change in terms of liquid fraction-temperature relationships. Charging and discharging packed bed column experiments have been also carried out for different operating conditions to analyze the potential of GPCC for packed bed thermal energy storage. The present column results clearly demonstrate the dependence of temperature variation along the packed bed and the overall performance of the storage unit on the phase change characteristics of GPCC. Small and non-uniform particles diameters of GPCC and heterogeneity of the bed material complicate the phenomena of heat transfer and evolution of latent heat in the packed bed. Mathematical modeling of the packed bed that considers the GPCC and air as two separate phases with inter-phase heat transfer is presented. Comparisons between experimental and numerical results are used to evaluate the sensitivity of numerical simulations to different model parameters.

Introduction

Research on advanced materials and systems for latent heat thermal energy storage (LHTES) has been intensified in recent years as due to the increasing interest in the design and development of energy-efficient building and the greater use of natural and renewable energy sources for heating, ventilation and cooling applications [START_REF] Farid | A review on phase change energy storage: materials and applications[END_REF]. Different types of storage units have been conceived and studied. They differ by the nature of phase change material (PCM) and geometry of the unit. Storage of thermal latent heat in packed bed of spheres, containing each a PCM, provides direct heat exchange between the PCM and the heat transfer fluid. The analysis of packed beds thermal energy storage has been the subject of intensive research work [START_REF] Ismail | Numerical and experimental study of spherical capsules packed bed latent heat storage system[END_REF][START_REF] Benmansour | Experimental and numerical investigation of solid particles thermal energy storage unit[END_REF]. A recent review on the heat transfer characteristics of thermal energy storage systems using PCM capsules [START_REF] Regin | Heat transfer characteristics of thermal energy storage system using PCM capsules: A review[END_REF] shows that most of the previous studies have been limited to packed beds of spheres or cylinders containing PCM with large characteristic diameter in the order of few centimetres. In recent years, granular phase changing composites (GPCC) obtained by micro or macro-encapsulating PCMs in highly porous solid structures with protecting envelopes have been developed using different encapsulation techniques [START_REF] Zhang | Granular phase changing composites for thermal energy storage[END_REF]. A wide range of GPCC in the form of particles of small diameters in the order of few millimetres are available that are suitable for low temperature energy storage applications.

Utilization of granular phase change composites (GPCC) of small particle diameter (1-3 mm) in latent heat thermal energy storage (LHTES) systems is the main concern of the present study. The motivations of the present interest in these materials are related to expected advantages in terms of heat transfer rates, flexibility, and multiplicity of applications in energy systems offered by using such GPCCs of small particle diameter (see Section 7 for a general overview). Potentials of their application have been recently investigated in desiccant cooling systems [START_REF] Rady | Study of heat and mass transfer in a dehumidifying desiccant bed with macro-encapsulated phase change materials[END_REF] and in air conditioning systems of buildings [START_REF] Nagano | Study of a floor supply air conditioning system using granular phase change material to augment building mass thermal storage-Heat response in small scale experiments[END_REF].

RUBITHERM GR27 and GR41 GPCC materials have been chosen for the present investigation [8]. They are considered as two examples representative of GPCC materials and the results obtained in the present study are useful in evaluating the application of such materials in LHTES systems and highlighting difficulties associated with characterization and modeling. The composition of the granule is 65% ceramics and 35% paraffin wax by weight with 1-3 mm particle diameter. GPCC materials offer the advantage of maintaining their macroscopic solid form during phase change. The PCM is bound within a secondary supporting structure (polymer or ceramic) that ensures that the PCM, when in the liquid form, does not leak out of the granulate. The result is that the bound PCM is always a "dry" solid product and liquid handling is eliminated.

In the process development stage of a LHTES system, the knowledge of thermophysical properties of materials used for latent heat storage is essential for design and optimization. The The T-history method is an alternative simple and low cost technique for characterization of PCMs. In the original T-history method [START_REF] Zhang | A simple method, the T-history method, of determining the heat of fusion, specific heat and thermal conductivity of phase-change materials[END_REF], the PCM sample and a sample with known thermal properties (distilled water) of the same initial homogeneous temperature are suddenly subject to ambient air with a room temperature which can be time dependent. Their temperature history upon cooling down is recorded.

Comparison of both curves, using a mathematical description of the heat transfer, allows the determination of the heat capacity of the solid and liquid phases and the latent heat of fusion of the PCM sample from the known specific heat of the reference material.

Using the T-history method, further modifications and improvements are employed to handle granular materials undergoing phase change over a temperature range. The accuracy of the original T-history method is shown to be limited by the assumption of constant temperature independent specific heats and the difficulty of determining the limits of solid and liquid phases. The concept of enthalpy and its relationship with temperature has been employed in the analysis to overcome these difficulties. Enthalpy-temperature and apparent heat capacity curves similar to those obtained using DSC have been developed for the GPCC, see Figs. 1 and2. Characterization results for GR27 and GR41 are shown in Table 1. These include the start and end temperatures of phase change (T m1 , T m2 ), the specific heat of solid and liquid phases (c ps , c pl ), and the latent heat of phase change (H ls ) which is the difference in enthalpy at temperature T m2 and T m1 . The values of T m1 and T m2 are obtained from the apparent heat capacity curve in a manner similar to the analysis of the heat fluxtemperature curve using DSC. assumed to be negligible. The air and GPCC energy equations are written as [START_REF] Beasley | Transient response of a packed for thermal energy storage bed[END_REF]:

( ) ( ) a a pa a a d pa a a a a sa sa s a c T U c T T k h S T T t z z z ∂ ρ ε ∂ρ ∂ ∂ + = ε + - ∂ ∂ ∂ ∂ (1)
c T T f k h S T T H t z z t ∂ ρ ε ∂ ∂ ∂ = ε - - -ρ ε ∂ ∂ ∂ ∂ (2)
The inter-phase heat transfer is described by the heat transfer coefficient (h sa ). The evolution of latent heat during phase change of the GPCC is accounted for by using the source

term s ls s l H f t ρ ε ∂ ∂ in Eq. ( 2

ACCEPTED MANUSCRIPT

aeff a d k k k = + (5)
The dispersive component of the longitudinal thermal conductivity in the porous media was evaluated as function of the Peclet number [START_REF] Kuwahara | A numerical study of thermal dispersion in porous media[END_REF][START_REF] Amiri | Transient analysis of incompressible flow through a packed bed[END_REF]

= < -ε = > ε (6)
To complete the mathematical description of the present problem, the initial and boundary conditions are specified as follows: at t=0: T a =T s =T o , at z=0:

T a =T in , = ∂ ∂ z T s / 0, at z=L: / / s a T z T z 0 ∂ ∂ = ∂ ∂ = .
Where T o is the initial bed temperature and T in is the inlet air temperature.

Numerical Solution

Examination of the governing equations ( 1) and ( 2) shows that the liquid fraction f l represents an extra unknown that is known only after the solution for solid temperature T s . Therefore, an iterative technique has been adopted for the calculation of f l as follows, ( )

ps s m i 1 i l l f c T T f f where the correction H + - = + ξ ξ = ( 7
)
Where i is the iteration number and the constraints max . , The value of T m used in Eq.( 7) is therefore function of the liquid fraction T m =f(f l ). This functional relationship is a characteristic of the phase change material as explained afterwards.

Equations ( 1) and ( 2) are discretized by using the finite volume method employing an implicit time scheme. The upwind scheme is used for discretization of the convection terms and the central difference scheme is used for discretization of the diffusion terms. The resulting equations are solved using the Gauss-Seidel iteration method. A time step of 0.001 s and a grid mesh of 0.1 mm are used in the simulations to ensure numerical stability and accuracy. The simulation results are independent of the grid size and time step used.

Figures 6 and7 show complete comparisons of experimental (T1, T8, T3, T9, T7, Tout) and numerical simulation (T1 Sim, T3 Sim, T7 Sim, Tout Sim) results during charging of GR27 and GR41 columns. In these calculations, Ranz model for particle-to fluid heat transfer has been used, axial dispersion is calculated using Eqs. ( 5) and ( 6), an average bed porosity of 0.42 is used, and the phase changing characteristics resulting from the T-history method are employed for the functional relationship of liquid fraction and Temperature. The agreement between experimental and numerical results is reasonable. However, some differences are observed and need to be evaluated.

These differences may be attributed due to uncertainty in the parameters of heat transfer model such as the phase change characteristics of GPCC, bed porosity, fluid-particle heat transfer coefficient, particle diameter, and axial dispersion. An analysis of the sensitivity of agreement between numerical and experimental results to these parameters is performed in next section in order to define the most critical model components.

Sensitivity of Agreement to Model Parameters

The sensitivity of numerical results to the values of particle-to-fluid heat transfer coefficient, dispersive component of thermal conductivity, and description of the phase change characteristics of GPCC are investigated. Numerical simulations have been carried out to evaluate the effects of these model components on the agreement between the experimental and numerical results. The effect of each parameter is investigated separately and unless otherwise specified the parameters used in the above simulations of Figures 6 and7 are used in the following simulations.

Effects of phase change characteristics:

The evolution of latent heat during charging and discharging of LHTES unit as function of temperature is an important issue for modeling, optimization, and design of the unit. This evolution is expressed in terms of the variation of liquid fraction f l of GPCC with temperature as follows.

( ) ( ) ( )

s m l ls H T H T f T H - = (8)
H(T) is the enthalpy at any temperature (T) and H s (T) is the solid phase enthalpy. Variations of f l with temperature for GR27 and GR41 are obtained by further analysis of T-history and DSC measurements. Applications of different forms for the functional variation of liquid fraction with temperature in numerical simulations of packed beds may affect the bed performance characteristics.

Applications of different functions of f l -T results in different evolution of latent heat during phase change. Simulations have been also carried out using an assumed linear profile of f l -T between T m1 and T m2 , an assumption that is usually used to handle phase change over a temperature range, and comparisons have been made with the results obtained using measured f l -T curves obtained from the T-history and DSC. Figure 8 shows that the assumption of linear profiles of f l -T can induce remarkable differences in the predicted variation of bed temperature with time. The differences can be clearly distinguished during the phase change period. The results obtained using curves of T-history and DSC are close to the experimental measurements.

Effects of bed porosity:

Packed bed porosity has been the subject of extensive research [START_REF] Theuerkauf | Analysis of particle porosity distribution in fixed beds using the discrete element method[END_REF][START_REF] De Klerk | Voidage variation in packed beds at small column to particle diameter ratio[END_REF]. Some of this research work has focused on developing correlations to predict overall void fraction for use in global prediction of pressure drop as function of the ratio of tube diameter to some effective particle diameter, D/d. The average bed porosity is modeled by deKlerk [START_REF] De Klerk | Voidage variation in packed beds at small column to particle diameter ratio[END_REF] as an exponential decay function.

. exp . using different values of bed porosity show that the effect of bed porosity is significant. The local changes in porosity can lead to large variations in predicted velocity profile and therefore nonuniform head loss along the packed bed. The accurate prediction of local porosity is therefore important for predicting heat transfer in packed beds. This is especially true for packed beds having a particle size and shape distribution as the present application. It seems that some form of multidimensional statistical description of the packing will be necessary for performance modeling.

The average bed to particle diameter employed in the present study 25 provides a central region essentially free from the wall effect on void fraction. However, near the walls, the void fraction increases and may result in flow channeling.

Effects of fluid-to-particle heat transfer coefficient and axial dispersion:

Uncertainty of the fluid-particle heat transfer coefficient is related to the correlations existing in the literature. For example, the widely used correlations given by Eqs.( 3) and ( 4) for Ranz model and Galloway and Sage [START_REF] Galloway | A model of the mechanism of transport in packed distended and fluidised beds[END_REF] result in heat transfer coefficient that differ by about 15%. Also, the importance of axial dispersion should be evaluated. Another source of uncertainty can be attributed as due to the assumption of uniform particle diameter. Screen analysis of GR27 and GR41 shows that they contain 0,1% of particles with diameter d p > 4 mm, 7,5% with d p = 3,0 -4,0 mm, 89,5%

with d p = 1,0 -3,0 mm, 2,3% with d p = 0,5 -1,0 mm, and 0,6% with d p < 0,5 mm. Therefore, a mean particle diameter of 2 mm has been assumed for the purpose of numerical simulations. The particle non-uniform shape may affect the boundary layer and probably increase the heat transfer coefficient.

Simulations have been carried out to evaluate the effects of uncertainty of these model parameters on the agreement between numerical and experimental results. It has been found that the utilization of Ranz model and Galloway and Sage [START_REF] Galloway | A model of the mechanism of transport in packed distended and fluidised beds[END_REF] correlation for estimation of the fluid-toparticle heat transfer coefficient has no remarkable effect on the agreement between numerical and experimental results. Also, simulations results using different particle diameter of 1.5 mm are very similar to those obtained using a particle diameter of 2 mm as an average value. This can be explained by the fact that decreasing particle diameter decreases the heat transfer coefficient as due to the decrease of Reynolds number. On the other hand, the surface area for heat exchange between the particles and fluid increases with the decrease of the particle diameter. The final effect, within the uncertainty limits of particle diameter, is a small variation of the heat exchange rates between the particles and fluid that is manifested in unremarkable differences in the agreement between the experimental and numerical results. The effects of axial dispersion have been also evaluated by comparison of results obtained with and without axial dispersion. The results indicate that, for the range of Reynolds number studied, the longitudinal thermal dispersion has no significant contribution to the overall temporal thermal response of the packed bed.

General overview for LHTES applications

An analysis of the role of the characteristic diameter of PCM capsule in a packed bed thermal energy storage system is an essential issue that explains the importance, limitations, and advantages of utilization of GPCC of small characteristic diameter. The goal of an energy storage system is to obtain high heat transfer rates during charging and discharging with reasonable pumping requirements. Applying heat exchangers principles [START_REF] Incropera | Fundamentals of Heat and Mass Transfer[END_REF] and introducing the column diameter (D) as a characteristic length, the number of transfer units (NTU) can be written as:

6 (1 ) 6 (1 ) a a d a hA h D D NTU St m cp d U cp d ε ε ρ - = = = - ( 10
)
Where m is the mass flow rate of air, A is the total surface area of heat transfer, h is the particle to fluid heat transfer coefficient, and St is the Stanton number (St=Nu/Re.Pr). The functional dependence of heat transfer coefficient on Reynolds number based on particle diameter (Re) is given by Eq. ( 4). NTU is a dimensionless parameter whose magnitude influences the unit performance. High values of NTU result in enhancing the unit effectiveness with high charging and discharging rates. The pressure gradient along the bed length can be obtained from a modified form of the Ergun equation [START_REF] Ergun | Fluid flow through packed columns[END_REF][START_REF] Macdonald | Flow through porous mediathe Ergun equation revisited[END_REF].

d d p U U z d d µ ε ρ ε ε ε ∂ - - - = + ∂ (11)
Equations ( 11) can be written in dimensionless form by introducing the definition of friction factor (f) and considering the pressure drop per unit length of column ( P/L) as: 2 3

( / ) (1 ) 150(1 ) 1.75 Re

d P D L D f U d ε ε ρ ε ∆ - - = = + (12)
Equations ( 10) and ( 12 (10) to estimate NTU and Eq. ( 12) to calculate P(D/L). Using GPCCs of small particle diameter may be useful when limited volumes of energy storage are applicable. The storage granulate offers the clear advantage that it can be used to fill containers of any conceivable geometry (flexibility in applications). The granulate form of the latent heat storage material provides the large surface area required for effective heat transfer. Mathematical modeling of the packed bed that considers the GPCC and air as two separate phases with inter-phase heat transfer is presented. Comparisons between experimental and numerical results

CONCLUSION

show an important sensitivity of their agreement to the phase change characteristics of GPCC and estimation of bed porosity. The sensitivity of agreement to the estimation of particle-to-fluid heat transfer coefficient and the axial dispersion are found to be less important. The functional variation of liquid fraction of GPCC with temperature during phase change is therefore an important parameter that controls the dynamics of packed bed during charging and discharging processes. Also, for packed beds having a particle size and shape distribution as the present application, accurate prediction of local porosity is essential for analysis of heat transfer in the bed. It seems that some form of multidimensional statistical description of the packing will be necessary for performance modeling.

The present analysis is useful for practical applications in the design and optimization of LHTES system using GPCC during the process development stage and for performance evaluation of packed bed latent heat storage columns. In general, utilization of encapsulated PCM of small particle diameter is advantageous in terms of heat transfer rates and pumping requirements for applications involving low values of Reynolds number, compact systems, and systems with limited volume of energy storage.

List of Figures and Tables:

Table 1. Characterization of GR27 and GR41 using the T-history method.

Fig. 1 Enthalpy and apparent specific heat variation with temperature for GR27.

Fig. 2 Enthalpy and apparent specific heat variation with temperature for GR41.

ACCEPTED MANUSCRIPT

d

  Particle diameter (m) T m phase change temperature (K) D packed bed column diameter (m) U d Darcy velocity (m/s) f friction factor Z distance along the bed height (m) f l liquid fraction Greek symbols G γ mass flux [kg/(m 2 s)] ε porosity; volume fraction h sa fluid-to-particle heat transfer coefficient [W/(m 2 K)] µ air dynamic viscosity (kg/m s) H enthalpy (J/kg) ν air kinematic viscosity (m 2 /s) H ls PCM latent heat of fusion (J/kg)

  quantity of thermal energy possible to store depends on the enthalpy variation in the working temperature range. Also, the enthalpy function or the apparent heat capacity of the PCM has to be known in order to obtain sufficiently accurate results in numerical simulations. This should take into account the temperature range over which the phase transition occurs under real operating conditions. Difficulties in the characterization of GPCC are related to the specific nature of the material and heterogeneity of the sample and require special care and improvements in the analysis of measurement results. Moreover, small and non-uniform particles diameters of GPCC and heterogeneity of the bed material complicate the phenomena of fluid-particle heat transfer and evolution of latent heat in the storage packed bed. Studies on the utilization of GPCC taking into account these phenomena are important for successful integration in LHTES systems. The present study aims at providing experimental data on the characterization of phase changing parameters of GPCC and highlights the difficulties and required improvements in the measurement methods. Packed bed column experiments during charging and discharging modes have been carried out for different operating conditions to analyze the potential of GPCC for packed bed thermal energy storage and provide experimental data for comparison with numerical simulation. Comparisons between experimental and numerical results have been also carried out to show the sensitivity of their agreement to parameters concerning phase change characteristics of GPCC, bed porosity, heat transfer correlations between the particle and fluid, and axial dispersion.2. Materials CharacterizationDifferential scanning calorimetry (DSC) and a modified T-history method have been used to study the phase changing behavior of granular composites. Details on the experimental set-up and data analysis for GR27 and GR41 can be found in[START_REF] Rady | A comparative study of phase changing characteristics of granular phase changing materials using DSC and T-history methods[END_REF]. Characterizing granulate PCMs using DSC requires special care to assure that the thermophysical properties of samples (1-10 mg) are those of the bulk materials. The main difficulty of DSC characterization of granular PCM lies in the heterogeneity and structure of the sample as it significantly affects the details of the heat transfer between the sample pan and PCM material. Direct utilization of the measured DSC curves could result in an inexact representation of the sample enthalpy change. A simple procedure has been advised to obtain accurate results from the DSC measurements based on the estimation of the thermal resistance between the sample and its enclosure[START_REF] Rady | A comparative study of phase changing characteristics of granular phase changing materials using DSC and T-history methods[END_REF]. The thermal resistance between the GPCC sample and it container is evaluated using DSC measurements obtained at three different values of cooling/heating rates. Estimation of the thermal resistance allows correcting the measured DSC heat flux and calculation of corresponding correct sample temperature. The corrections are found to be significant at high values of heating/cooling rates and within the phase change temperature range of GPCC. High DSC heating or cooling rates promote the development of thermal and chemical non-equilibrium effects inside the sample material and result in an increase in the temperature range for melting and solidification. The latent heat of phase change is found to be essentially independent on the DSC heating/cooling rate.

Figures 4

 4 Figures 4 and 5 show the time variation of temperature of the inlet air and at 4 cm from the

  ). Where f l is the liquid fraction of GPCC undergoing phase change and H ls is the latent heat of fusion of the PCM material. The value of f l is related to the temperature s T and determines the rate of evolution of latent heat during phase change.Knowledge of fluid-to-particle heat transfer coefficient sa h is an important issue for the analysis of heat transfer in packed beds. A considerable amount of research work has been carried out in the literature to identify suitable correlations for the particle to fluid heat transfer coefficient.Numerous correlations for predicting heat transfer coefficients in packed beds exist in the literature[START_REF] Baptista | Dimensionless analysis of fluid-to-particle heat transfer coefficients[END_REF][START_REF] Sanderson | Performance and efficient design of packed bed thermal storage systems[END_REF][START_REF] Galloway | A model of the mechanism of transport in packed distended and fluidised beds[END_REF][START_REF] Wakao | Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds[END_REF][START_REF] Nosfor | Measurements of the gas-particle convective heat transfer coefficient in a packed bed for high temperature energy storage[END_REF][START_REF] Seshadri | Comparison of formulae for determining heat transfer coefficient of packed beds[END_REF]. A simple analysis of heat conduction from an isolated spherical particle surrounded by a stagnant fluid leads to the conclusion that the limiting value of Nusselt number for low Reynolds number flows is 2.0. At higher Reynolds convective effects become significant and a correlation of the form Nu=f(Pr, Re) is to be expected. Correlations based on the Ranz model[START_REF] Ranz | Friction and transfer coefficients for single particles and packed beds[END_REF] and Galloway and Sage[START_REF] Galloway | A model of the mechanism of transport in packed distended and fluidised beds[END_REF] have been widely used for heat transfer modelling in packed beds[START_REF] Benmansour | Experimental and numerical investigation of solid particles thermal energy storage unit[END_REF][START_REF] Beasley | Transient response of a packed for thermal energy storage bed[END_REF]. The Ranz model provides a method of applying convection correlations for single spheres in an external flow to spheres in packed beds. With the Ranz model, the superficial velocity U d is multiplied by a factor of 10.73 to yield an effective free stream velocity inside the packing. Thus the Reynolds number velocity of 10.73U d is used with convection correlations for single spheres to predict convection coefficients in packed peds. Ranz determined the factor 10.73 analytically and verified it by comparison with experimental results. The Nusselt number correlation using Ranz model [18] is written as: [19] developed the following correlation for fluid-to-particle heat transfer coefficient in packed beds. . For randomly packed bed of spheres c 1 and c 2 are 1.354 and 0.0326, respectively. Beasley and Clark [11] used higher values of these coefficients (c 1 =2.031 and c 2 =0.049 ), resulting in an increase of about 50% in the heat transfer coefficient, and obtained better agreement between the numerical and experimental results.In the application of the above expressions for heat transfer, it is important to note that the axial dispersion coefficients as due to axial mixing must be properly estimated. There are two main mechanisms which contribute to axial dispersion: molecular diffusion and turbulent mixing arising from the splitting and recombination of flows around the particles. The effective thermal conductivity of the fluid phase consists of the stagnant and the dispersion conductivity.

  At convergence =0 and the temperature of phase changing nodes is equal to the phase change temperature T m . It should be mentioned that the correction term takes the form of Stefan number. Previous simulations of packed beds reported in the literature have been concerned with phase changing materials at constant temperature T m . In such cases, the value of f l changes from zero to one during melting at T s =T m . However, the present characterization of GPCC shows that the phase change process takes place over a temperature range (between T m1 and T m2 ).

  the packing mode is embodied in the average bed porosity at infinite column diameter b . Literature values vary for the bed voidage of a randomly packed bed of spheres at large column to particle diameter ratio. This variance is ascribed to the packing mode. Four random packing modes for spheres are distinguished: (a) very loose random packing ( b 0.44) obtained by gradual defluidization of a fluidized bed or by sedimentation; (b) loose random packing ( b 0.40-0.41) obtained by letting spheres roll individually in place, or by dropping the spheres into the container as a loose mass; (c) poured random packing ( b 0.375-0.391) obtained by pouring spheres into a container; and (d) dense random packing ( b 0.359-0.375) obtained by vibrating or shaking down the packed bed. Due to the stochastic nature of the packing process, accurate and reliable predictions of bed porosity are elusive. The sensitivity of numerical predictions to the value of bed porosity is shown in Fig. 9. Comparisons of experimental and numerical results obtained

  ) are used to investigate the variation of NTU and pressure drop along the bed considering air as the heat transfer fluid at different values of PCM capsule diameter. In order to make this analysis, one has to consider a constant Reynolds number based on the column diameter different values of D/d, calculate (Re=Re o d/D), and then employ Eq.

Figure 10 shows

 10 Figure 10 shows the variation of NTU and pressure drop along the bed with D/d for different

  Utilization of granular phase change composites (GPCC) of small particle diameter(1-3 mm) in latent heat thermal energy storage (LHTES) systems is the main concern of the present study. Experimental measurements using DSC and T-history methods are carried out for characterization of phase changing parameters of GPCC such as phase change temperature, latent heat, and energy storage capacity. The evolution of latent heat with temperature during phase change in terms of liquid fraction-temperature relationships is described. Difficulties and required improvements in the measurement methods to handle GPCC are discussed. Packed bed column experiments during charging and discharging modes are carried out for different operating conditions to analyze the potential of GPCC for packed bed thermal energy storage and provide experimental data for comparison with numerical simulation.The present column results clearly demonstrate the dependence of temperature variation along the packed bed and the overall performance of the storage unit on the phase change characteristics of GPCC. Small and non-uniform particles diameters of GPCC and heterogeneity of the bed material complicate the phenomena of heat transfer and evolution of latent heat in the packed bed.

Fig. 3

 3 Fig. 3 Schematic of the experimental set-up for packed-bed.

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4 Variation of bed temperature during charging and discharging of GR27 column, Re =7.29.Fig. 5 Variation of bed temperature during charging and discharging of GR41 column, Re =7.29.Fig. 6 Comparison of experimental and numerical results, (GR27, Re=7.29).

Fig. 7

 7 Fig. 7 Comparison of experimental and numerical results, (GR41, Re=7.29).

Fig. 8

 8 Fig. 8 Comparison of experimental and numerical simulations using different liquid fraction functions, (GR27, Re=7.29).

Fig. 9

 9 Fig. 9 Comparison of experimental and numerical results, effects of bed porosity, (GR27, Re=7.29).

Fig. 10

 10 Fig. 10 Variation of NTU and pressure drop with PCM capsule diameter.

Figure 1 :

 1 Figure 1: Enthalpy and apparent specific heat variation with temperature for GR27.

Figure 2 :

 2 Figure 2: Enthalpy and apparent specific heat variation with temperature for GR41.

Figure 3 :

 3 Figure 3: Schematic of the experimental set-up for packed-bed.

Figure 4 :

 4 Figure 4: Variation of bed temperature during charging and discharging of GR27 column, Re =7.29.

Figure 5 :

 5 Figure 5: Variation of bed temperature during charging and discharging of GR41 column, Re =7.29.

Figure 6 :

 6 Figure 6: Comparison of experimental and numerical results, (GR27, Re=7.29).

Figure 7 :

 7 Figure 7: Comparison of experimental and numerical results, (GR41, Re=7.29).

Figure 8 :

 8 Figure 8: Comparison of and numerical simulations using different liquid fraction functions, (GR27, Re=7.29).

Table 1 :

 1 Characterization of GR27 and GR41 using the T-history method.

	Sample	T m1 (°C)	T m2 (°C) c ps (J/kg.K)	c pl (J/kg.K)	H ls (J/kg)
	GR27	21	29	2044	1921	64850
	GR41	31	45	2784	2745	65895


F Benatti
email: benatti@ts.infn.it

R Floreanini

U Marzolino

Entanglement and Squeezing with Identical Particles: Ultracold Atoms Quantum Metrology

Entanglement and squeezing with identical particles:

ultracold atom quantum metrology

Ultracold gases trapped in optical lattices provide among the most successful physical setups for studying quantum many-body effects in condensed matter physics [1] and quantum information [2,3]; more recently, these systems have been proposed as experimental realizations of quantum metrological protocols [4,5] In all these applications, the property of entanglement [6] plays a crucial role; no longer an epistemological curiosity, entanglement is indeed becoming more an more an experimentally accessible resource, with applications in many areas of many-body physics. In particular, entangled N -qubit states have been proposed as means to beat the so-called shot-noise limit accuracy in parameter estimation [7,8].

A first step in this direction using many-body systems has been recently realized: entangled states in systems of ultra-cold atoms have been generated through spinsqueezing techniques [9,10], with the aim of using them as input states in interferometric apparatuses, specifically constructed for quantum enhanced metrological applications. In such devices, the initial N -qubit states are rotated by means of collective spin operators. In the case of distinguishable qubits, all these rotations are local; therefore, preliminary spin squeezing of separable states is necessary in order to introduce (nonlocal) quantum correlations, thus opening the possibility of sub-shot-noise accuracies in parameter measurements.

For systems of distinguishable particles, the notion of quantum entanglement is well-established [6]: a generic state of N qubits is defined to be entangled if it is not fully separable, namely if it cannot be written as

ρ sep = k p k ρ (1) k ⊗ ρ (2) k • • • ⊗ ρ (N ) k , (1)
Confidential: not for distribution. Submitted to IOP Publishing for peer review 1 April 2011

where the p k ≥ 0 are weights, k p k = 1, while each ρ

(j)
k is a density matrix for the j-th qubit acting on the corresponding Hilbert space À j .

However, in the case of trapped ultra-cold atoms, the qubits involved are identical, a fact that has often not been fully appreciated in the recent literature dealing on quantum metrology. Indeed, the indistinguishability of the particles asks for a rethinking of definitions and concepts originally introduced in the case of distinguishable qubits. More specifically, the notion of separability based on Eq.( 1) is strictly associated with the tensor product structure of the Hilbert space, À N = N j=1 À j , which is natural in the case of N distinguishable particles. On the other hand, pure bosonic states must be symmetric under exchange of particles and mixed states must be convex combinations of projections onto such states. This fact demands a different approach to the notions of non-locality and entanglement based not on a structure related to the particle aspect of first quantization, as in (1), rather on the behavior of correlation functions of commuting observables more generally related to the mode description typical of second quantization. ‡ In this spirit, a generalized notion of entanglement has been introduced in [14] which reduces to the standard one for distinguishable qubits. Although the notion of entanglement in many-body systems has already been addressed and discussed in the literature (for instance, see [2], [11]- [13], [15]- [27] and references therein), only limited results actually apply to the case of identical particles; the new definition of separability in [14] is instead very general, suitable to all physical situations and notably to systems of trapped ultracold atoms.

In the following, we show that, according to such new definition, neither entangled states, nor spin-squeezing are necessary in order to achieve sub-shot noise accuracies in quantum metrological applications based on systems of identical particles. Though necessary, quantum non-locality comes not from the states, rather from the rotations that are implemented in the interferometric apparata. Indeed, while all these rotations are local for distinguishable qubits, in the case of identical qubits some of those operations, implemented through the action of collective spin operators, happen to be non-local and might allow sub-shot noise accuracies without the need of a preliminary squeezing.

In general, density matrices as in (1) are not allowed bosonic identical particle states, not even when the single-particle states are all equal in each term of the sum, i.e. ρ 2. Since Ψ -|ρ ⊗ ρ|Ψ -= Det(ρ), the density matrix ρ ⊗ ρ cannot correspond to a state of two bosonic qubits unless ρ is a projection (i.e. Det(ρ) = 0).

Therefore, the tensor product structure which is natural for distinguishable particles is not appropriate for discussing the entanglement properties of systems of identical particles. These should rather be investigated within the second quantization formalism whereby one introduces creation a † i and annihilation operators a i of single-particle orthonormal basis states |i , obeying the canonical commutation relations [a i , a † j ] = δ ij . Entanglement in such a context should correspond to whether, given a state ω of the system, there are non-classical correlations among commuting observables [14] as, for instance, those localized in non-overlapping spatial regions V 1 and V 2 [12]. Specifically, this can be inspected by considering the structure of two-point functions of the form ω(P V 1 P V 2 ), where P Vα , α = 1, 2, are polynomials in the creation and annihilation operators a † (ψ α ), a(ψ α ), of states |ψ α spatially localized in the region V α . The relevant fact is that the spatially local algebras generated by the polynomials

P Vα commute: [a(ψ 1 ) , a † (ψ 2 )] = ψ 1 |ψ 2 = 0.
More in general, one can argue about the entanglement between observables belonging to two generic commuting sub-algebras (A, B) of the entire algebra generated by creation and annihilation operators, which we shall refer to as algebraic bipartition [14]. We shall call an observable (operator) O "(A, B)-local" if it is of the form O = AB, A ∈ A and B ∈ B, and a state ω "(A, B)-separable" if the expectation values ω(AB) of local operators AB can be decomposed into a convex linear combination of product of expectations:

ω(AB) = i λ i ω a i (A) ω b i (B) , λ i > 0 , i λ i = 1 , (2)
in terms of other states ω a i , ω b i ; otherwise, ω is said (A, B)-entangled. Remark 1: Algebraic and standard entanglement. As already mentioned, the previous definition reduces to the standard notion of entanglement for distinguishable particles [14]. This becomes apparent in the case of two qubits, by choosing the algebraic bipartition A = B = M 2 , where M 2 is the algebra of 2 × 2 matrices over 2 , and the expectation value ω(AB) = Tr (ρ A ⊗ B) defined in the usual way through the trace with the two-qubit density matrix ρ.

Let us stress that in the case of identical particles, there is no a priori given bi-partition so that questions about entanglement and separability, non-locality and locality are meaningful only with reference to a specific class of (commuting) observables; the general definition given above precisely stems from this fundamental observation.

A concrete application of the previous considerations is the second quantization of a single-particle with Hilbert space 2 which, in the Bose-Hubbard approximation, effectively describes N ultracold atoms confined by a double-well potential. Then, the state | ↓ describes one atom located within the left well and the state | ↑ an atom localized within the right one. Let |0 be the vacuum state and a † , b † the creation operators of a particle in the states

| ↓ and | ↑ , that is a † |0 = | ↓ , b † |0 = | ↑ .
When the total number N is conserved, the symmetric Fock space of this twomode system is generated by N + 1 orthonormal eigenvectors of the number operator

N = a † a + b † b: |k = (a † ) k (b † ) N -k k!(N -k)! |0 , 0 ≤ k ≤ N . (3)
Because of the orthogonality of the spatial modes, by considering the norm-closures of all polynomials P a in a, a † , respectively P b in b, b † , one obtains two commuting subalgebras A and B.

According to the definition in 2), the states |k are (A, B)-separable; indeed, they are created by the (A, B)-local operators (a † ) k (b † ) N -k . More in general, one can show that (A, B)-separable density matrices must be convex combinations of projections |k k| [14]:

ρ = N k=0 p k |k k| , p k > 0 , N k=0 p k = 1 . (4)
Consider instead the following operators,

J x = 1 2 (a † b + ab † ) , J y = 1 2i (a † b -ab † ) , J z = 1 2 (a † a -b † b) , (5)
J x = 1 2 (c † c -d † d) , J y = 1 2i (d † c -dc † ) , J z = 1 2 (c † d + cd † ) .
Relatively to (C, D), it is now e iθ Jx = e iθ c † c e -iθ d † d which acts locally.

Remark 2: Mode entanglement vs. particle entanglement. In first quantization, an (A, B)-separable state for a N = 2 bosonic qubits like

a † b † |0 corresponds to (| ↑↓ + | ↓↑ )/ √ 2.
Such state is surely entangled for distinguishable qubits, while, according to our definition, it is no longer so for identical bosonic qubits; the reason is that its entanglement is only formal as it comes from the necessary symmetrization of the separable state | ↑↓ [15]. A Bogolubov transformation as the one above corresponds to a change of basis in the single particle Hilbert space, from the one of spatially localized states, to the one of

c † |0 = (| ↓ + | ↑ )/2, d † |0 = (| ↓ -| ↑ )/2.
Physically speaking, such states are eigenstates of the single particle Hamiltonian in the Bose-Hubbard approximation with a highly penetrable barrier. The change to the energy bipartition (C, D) is non-local with respect to the spatial bipartition (A, B), though it corresponds to a local unitary rotation in first quantization.

In the case of a system of N distinguishable qubits, the collective angular momentum operators J α , α = x, y, z and the corresponding rotations are sums of single particle spin operators,

J (i) α , α = x, , y, z, i. e. J α = N i=1 J (i)
α . These operators are local with respect to the tensor product structure in (1).

Based on this, the variance ∆ 2 J n of the collective spin J n = n x J x + n y J y + n z J z along the unit spatial direction n = (n x , n y , n z ), with respect to separable vector states

|Ψ = N i=1 |ψ i results [28, 29] ∆ 2 J n = Ψ|J 2 n |Ψ -Ψ|J n |Ψ 2 = N 4 - N j=1 ψ i |J (i) n |ψ i 2 ≤ N 4 . (6)
Therefore, ∆ 2 J n is an entanglement witness for pure states, in the sense that if ∆ 2 J n > N/4 then the pure state |Ψ cannot be fully separable. This is no longer the case for N identical bosonic qubits. Indeed, consider the number states |k in (3); using (5), one gets

k|J n |k = n z 2 (2k -N ) ( 7
)
k|J 2 n |k = N + 2k(N -k) 4 + n 2 z N (N -1) -6k(N -k) 4 (8)
∆ 2 J n = 1 -n 2 z 4 N + 2k(N -k) . (9)
Therefore, if k = 0, N , for all n that satisfy

n 2 z < 2k(N -k) N + 2k(N -k) ≤ 1 ,
the states |k , though (A, B)-separable, nevertheless yield ∆ 2 J n > N/4; therefore, ∆ 2 J n is not an entanglement witness for pure states of bosonic qubits.

In greater generality, inequalities for mean values and variances of collective spin operators with respect to any (mixed) separable state of distinguishable qubits as in (1) have been derived in [30]; these are called spin squeezing inequalities and read: §

J 2 n 1 + J 2 n 2 + J 2 n 3 - N (N + 2) 4 ≤ 0 , (10)
∆ 2 J n 1 + ∆ 2 J n 2 + ∆ 2 J n 3 - N 2 0 , (11)
J 2 n 1 + J 2 n 2 - N 2 -(N -1)∆ 2 J n 3 ≤ 0 , (12)
(N -1)(∆ 2 J n 1 + ∆ 2 J n 2 ) -J 2 n 3 - N (N -2) 4 0 , (13)
where n 1 , n 2 , n 3 denotes any triplet of unit vectors corresponding to orthogonal spatial directions and X the mean value of an operator X.

It is thus interesting to study whether these inequalities are also satisfied by (A, B)-separable states (4) of N identical qubits. Denote by k = N k=0 p k k and k 2 = N k=0 p k k 2 the first and second moments of the N + 1-valued stochastic variable § In [30], these inequalities are derived with respect to the standard triplet n 1 = x, n 2 = ŷ, n 3 = ẑ. The result easily extends to more general triplets.

k with respect to the probability distribution π = {p k } N k=0 . Using ( 7) and ( 9), meanvalues and variances of the collective spin operators J n with respect to the states in (4) read

J n = n z 2 2 k -N (14)
J 2 n = N (1 + 2 k ) -2 k 2 4 + n 2 z 4 N (N -1) -6N k + 6 k 2 (15)
∆ 2 J n = N (1 + 2 k ) -2 k 2 4 + n 2 z 4 6 k 2 -2 k (N + 2 k ) -N . ( 16
)
From the orthogonality of the triplet n 1 , n 2 , n 3 , it follows that n 2 1z + n 2 2z + n 2 3z = 1; one can thus check that all the previous inequalities but (12) are satisfied by (A, B)-separable states. Concerning (12), its left hand side reads

δ = N 2 ∆ 2 k -k (N -k ) + n 2 3z 2 (N + 2) k (N -k ) -3N ∆ 2 k , (17)
where

∆ 2 k := k 2 -k 2 is the variance of k with respect to π = {p k } N k=0 . If π is chosen such that a := k (N -k ) > N∆ 2 k
, then δ becomes positive and thus (12) results violated by the corresponding (A, B)-separable states (4) for all orthogonal triplets with

N (a -∆ 2 k) (N + 2) a -3 N ∆ 2 k < n 2 3z ≤ 1 . (18)
Consider the pure states | in (3) with = 0, N ; in such a case, p k = δ k , k = and ∆ 2 k = 0, so that ( 12) is violated for N/(N + 2) < n 2 3z ≤ 1. Remark 3: Additional squeezing inequalities failing for identical qubits. Different spin-squeezing inequalities have also been presented in [31]. Although derived for distinguishable qubits, they have been directly applied to signal entanglement also in the case of systems of identical particles. A simple instance of their failure for (A, B)separable states is provided by inequality (10) in [31]: N ∆ 2 J n + J n 2 < N/4; indeed, it is always satisfied by any Fock state (3) with 0 < k < N when n is chosen along the z direction, thus erroneously signaling entanglement in such states.

The preceding results indicate that spin-squeezing inequalities that are derived for distinguishable qubits can not directly be used as entanglement witnesses in the context of identical qubits. Since the use of spin-squeezed states for metrological purposes have recently become the focus of much theoretical investigations [7,8,29,32], we now discuss the impact of particle indistinguishability on such an issue.

For any orthogonal triplet of space-directions n 1 , n 2 , n 3 , the Heisenberg uncertainty relations for the SU(2) operators J n read

∆ 2 J n 1 ∆ 2 J n 2 1 4 J n 3 2 . ( 19
)
One speaks of spin-squeezing when one of the variances can be made smaller than 1 2 J n 3 . The relevance of states satisfying this condition for achieving otherwise unavailable accuracies has been studied in relation to the measure of an angle θ by interferometric techniques. These are based on a rotation of an input state ρ into the final state

ρ θ = e -iθJ n 1 ρ e iθJ n 1 , (20)
and upon measuring on ρ θ the collective spin J n 2 , where n 2 ⊥ n 1 . By choosing the remaining orthogonal unit vector n 3 such that J n 3 = Tr(ρ J n 3 ) = 0, by error propagation, the uncertainty δθ in the determination of θ can be estimated by [8]:

δ 2 θ = ∆ 2 J n 1 ∂ θ J n 2 θ θ=0 2 = ∆ 2 J n 1 J n 3 2 = ξ 2 W N , ( 21
)
in terms of the spin-squeezing parameter

ξ 2 W := N ∆ 2 J n 1 J n 3 2 . ( 22
)
The value δ 2 θ = 1/N is called shot-noise limit; in the case of distinguishable qubits, it gives the lower bound to the attainable accuracies when the input state ρ is separable. Indeed, in such a case one finds ξ 2 W ≥ 1. This result follows from the inequality ξ 2 W ≥ ξ 2

S

where the new spin-squeezing parameter

ξ 2 S := N ∆ 2 J n 1 J n 2 2 + J n 3 2 , (23)
has been introduced in [32]; by means of the local structure of the collective spin operators J α , α = x, y, z, one can prove that ξ 2 S is always ≥ 1 for separable states of distinguishable qubits. Therefore, using distinguishable qubits, the shot-noise limit can be beaten, namely accuracies better than 1/N can be achieved only if ξ 2 S < 1, that is only by means of entangled states.

Let us instead consider N identical bosonic qubits in the (A, B)-separable pure states |k and any triplet of orthogonal spatial directions n 1 , n 2 , n 3 , with n 1 = ẑ. Using ( 7)-( 9) and n 2 1z + n 2 2z + n 2 3z = 1, one computes

ξ 2 S = N ∆ 2 J n 1 J n 2 2 + J n 3 2 = N 1 -n 2 1z n 2 2z + n 2 3z N + 2k(N -k) (2k -N ) 2 = N (N + 2k(N -k)) (2k -N ) 2 ≥ 1 , 0 ≤ k ≤ N . ( 24
)
In the case of the (A, B)-separable density matrices (4), first observe that, thanks to the Cauchy-Schwartz inequality, one has

N k=0 p k J n 2 k ≥ J n 2
, where X 2 k denotes the mean-value of X with respect to the number state |k . Then,

ξ 2 S ≥ N N k=0 p k ∆ 2 k J n 1 J n 2 2 + J n 3 2 = 1 -n 2 1z n 2 2z + n 2 3z N k=0 p k N N + 2k(N -k) N k=0 p k (N -2k) 2 ≥ N k=0 p k k -N 2 2 N k=0 p k k -N 2 2 ,
where the last inequality follow from the second line of (24). A further application of the Cauchy-Schwartz inequality to the right hand side of the last inequality yields ξ 2 S ≥ 1 for all (A, B)-separable states when n 1 = ẑ.

If one chooses n 1 = ẑ, in the case of (A, B)-separable mixed states, one finds ∆ 2 J z = 0 and J n 2 = J n 3 = 0; therefore, ξ 2 S (and ξ 2 W ) diverges. Instead, for (A, B)separable pure states |k also ∆ k J z = 0 whence ξ 2 S (and ξ 2 W ) are not defined and must thus be computed by means of suitable limiting procedures.

Let us consider the (A, B)-entangled vector state |Ψ = N k=0 √ p k |k , with real coefficients from a probability distribution π = {p k } N k=0 over the stochastic variable k. Then, from ( 7)-( 9) it follows that J y = 0 and ∆ 2 J z = ∆ 2 k; therefore ξ 2 S = ξ 2 W and

ξ 2 W = N ∆ 2 J z J x 2 = N ∆ 2 k N k=1 k(N -k + 1) √ p k p k-1 2 . (25)
In the case of a Gaussian distribution peaked around k = = 0, N ,

p k = 1 Z exp - (k -) 2 σ 2 , Z = N k=0 exp - (k -) 2 σ 2 , (26)
one finds

ξ 2 W = 2N + O e -1 2σ 2 ( + 1)(N -) + (N -+ 1) 2 . (27)
Thus, for sufficiently small σ, ξ 2 W < 1 for all = 0, N . On the other hand, by choosing

|Ψ = p N N k=0 , k =N/2 |k + (1 -p) |N/2 , 0 < p < 1 , (28)
it turns out that ∆ 2 k = p(N + 2)(N + 1)/12 so that (25) yields

ξ 2 W = N (N + 1) 12( √ 1 -p + q) 2 , where q = p N 2 (N + 2) k =N/2,N/2+1 k(N -k + 1) . Letting p → 0 one gets |ψ → |N/2 and, if N > 3, ξ 2 W → N (N + 1)/12 > 1 .
Remark 4: Drawbacks of spin-squeezing for identical qubits. The above two examples show that, when n 1 = ẑ, the spin-squeezing parameters are not well-defined: different values for ξ 2 W and ξ 2 S can be obtained by approaching a state |k via different limit procedures. This fact is also of practical importance: indeed, in [5] approximations to Fock states |k have been experimentally constructed that are characterized by spin-squeezing parameters ξ 2 W < 1. This property arises from the fact that those approximations are (A, B)-entangled states. The previous discussion shows that some care has to be taken in constructing the perturbations of |k ; indeed, not all (A, B)entangled states arbitrarily close to it automatically have ξ 2 W,S ≥ 1. Therefore, the spin-squeezing parameters ξ 2 W,S are not always useful for metrological applications, a better quantity is the so-called quantum Fisher information [33], which, as we shall show below, is continuous and well defined for all bosonic qubits.

In a measurement of the angle θ based on the state rotation (20), the error ∆θ given by a locally unbiased estimator of the angle θ is bounded by (see the Appendix)

∆ 2 θ ≥ 1 F [ρ, J n 1 ] , (29)
where F [ρ, J n 1 ] is the so-called quantum Fisher information associated with the rotation of ρ around n 1 .

In order to overcome the shot-noise limit ∆ 2 θ = 1/N , the quantum Fisher information must then be strictly larger than N . In the appendix it is also showed that, in full generality,

F [ρ, J n 1 ] ∆ 2 J n 2 ≥ J n 3 2 , (30)
where n 1 , n 2 , n 3 is a triplet of orthogonal spatial directions. Thus, if J n 3 = 0, one gets the following relation between the quantum Fisher information and the squeezing parameter ξ 2 W in ( 22):

F -1 [ρ, J n 1 ] ≤ ∆ 2 J n 2 J n 3 2 = ξ 2 W N . (31)
In the case of distinguishable qubits, from (29) and (31) it follows that spin-squeezing, namely ξ 2 W < 1, opens the possibility of achieving ∆ 2 θ < 1/N , thus of beating the shot-noise limit.

In the case of identical qubits and of (A, B)-separable states, the right hand side of the above inequality diverges when n 1 = ẑ, since J n 3 = 0, while it does not make sense if n 2 = ẑ for then also ∆ 2 J n 2 = 0; whence, as already observed, ξ 2 W is not defined. However, the quantum Fisher information is always well-defined. Indeed, using (37) in the Appendix and ( 9), one finds that, if k = 0, N ,

F [|k k|, J n ] = 4 ∆ 2 J z = (1 -n 2 z ) N + 2k(N -k) > N (32) for 0 ≤ n 2 z < 2k(N -k) N + 2k(N -k) < 1 . (33)
In particular,

F [|k k|, J y ] = N + 2k(N -k) > N ,
for k = 0, N , so that, according to (29), the (A, B)-separable pure states |k might overcome the shot-noise limit; the largest value is obtained with the twin Fock state |N/2 yielding F [ρ, J y ] = O(N 2 ), thus allowing to approach the so-called Heisenberg limit ∆ 2 θ = 1/N 2 .

Remark 5: Approximations of Fock states. Notice that even approximating a number state | by an experimentally more realistic superposition | , σ of states |k with coefficients as in (26), may beat the shot noise limit. Indeed, one computes

F [| , σ , σ|, J y ] = 4 ∆ 2 J y = N + 2 (N -) + O e -1 σ 2
, which can be kept > N by choosing a suitably small σ.

Instead, making the quantum Fisher information larger than N is impossible without (A, B)-non-locality; indeed, if both the state and the rotation operation are (A, B)-separable, one finds:

F [| |, J z ] = 4 ∆ 2 J z = 0. Even considering the (A, B)- entangled perturbation | , σ does not help, since F [| , σ , σ|, J z ] = 4 ∆ 2 J z = 8e -1 σ 2 + O e -2 σ 2
, making the lower bound to the error in (29) arbitrarily large as σ → 0.

When dealing with (A, B)-separable mixed states (4), by means of the equation (36) of the Appendix, one computes [14]:

F [ρ, J n ] = (1 -n 2 z ) N + 2N k -k 2 -4 N k=0 p k p k+1 p k + p k+1 (k + 1)(N -k) . (34)
Thus, if (32) holds for a certain | , then, by continuity, F [ρ, J n ] > N for any probability distribution π = {p k } N k=0 suitably peaked around k = . Hence, by feeding the interferometer appartus even with a mixed separable state on might still be able to overcome the shot-noise limit. On the other hand, using the squeezing parameter criteria, we would wrongly discard such states as not useful for metrological applications, since , as shown before, we know that for all such mixed states one has: ξ 2 W ≥ ξ 2 S ≥ 1. Conclusions. In the case of metrological applications based on ultracold atom systems, entangled states are not necessary to reach sub-shot noise accuracies in parameter estimation; indeed, as shown above, one can feed quantum interferometers with separable states of N identical qubits and yet obtain a value greater than N for the associated quantum Fisher information. This phenomenon is surely due to the non-local character of the system; however, the non-locality is not in the input states, which in fact are (A, B)-separable, rather in the (A, B)-non-local character of the staterotation (20) generated by J n = J z in the experimental apparatus. Were the qubits distinguishable, neither the state nor the rotation would carry elements of non-locality so that in order to beat the shot-noise limit, the state should be turned into an entangled one before sending it to input port of the interferometer. This is exactly what is done via a spin-squeezing technique in the experiments reported in [9,10].

Instead, the main point we make here is that in experiments involving identical qubits, no state entanglement, as for instance that obtained through a preliminary squeezing operation, is needed in order to beat the classical shot-noise accuracies. Even separable Fock states as in (3) could be used, the necessary quantum non-locality being provided by the measuring procedure itself via the mode rotation operated by the interferometric device. In line of principle, a twin-Fock state |N/2 can approach the Heisenberg limit in phase estimation. Furthermore, although the experimental preparation of Fock states may be problematic, we have also shown that sub-shotnoise accuracies can be reached by states sufficiently close to Fock states as those obtained in [5]. Therefore, the results presented here open new perspectives in the actual realization of cold atom based interferometers able to outperform any apparatus so far constructed for measuring ultra-weak effects.

Appendix

A most used definition of quantum Fisher information F [ρ, J n ] is given by [34] F [ρ,

From such an expression one sees that the quantum Fisher information is a continuous function of the state ρ and that, for pure states, Analogously, inequality (30) follows from the fact that Tr(ρ θ (J n 2 -J n 2 θ ) = 0 implies Tr ∂ θ ρ θ (J n 2 -J n 2 θ = ∂ θ J n 2 θ = J n 3 θ .

=

  ρ k for all i. Indeed, consider two qubits and fix an orthonormal basis C 4 |ij , i, j =↑, ↓. If they are bosons, their states should have vanishing components on the anti-symmetric state |Ψ -= (| ↓↑ -| ↑↓ )/ √

that satisfy the SU( 2 )

 2 algebraic relations [J x , J y ] = i J z and their cyclic permutations. They are all non-local with respect to the algebraic bipartition (A, B) and such are the exponentials e iθJx and e iθJy , while e iθJz = e iθa † a e -iθb † b is (A, B)-local. By means of a Bogolubov transformation to other creation and annihilation operators (c † , c) and (d † , d), such that a = (c + d)/ √ 2 and b = (c -d)/ √ 2, together with their hermitian conjugate, one obtains another bipartition (C, D) and rewrites

  J n ] := tr ρ L 2 , (35)where L, known as symmetric logarithmic derivative,∂ θ ρ θ θ=0 = (ρ L + L ρ)/2 = -i [J n , ρ]. Given a spectral decomposition ρ = j r j |r j r j |, one computes F [ρ, J n ] = 2 i,j : r i =r j (r i -r j ) 2 r i + r j r i |J n |r j 2 .

F 2 ≤

 2 [|ψ ψ|, J n ] = 4 ∆ 2 ψ J n . (37) An estimator E is locally unbiased if ∂ θ Tr(ρ θ E) θ=0 = 1; then, inequality (29) follows from applying to this relation the Cauchy-Schwartz inequality for matrices Tr(AB) Tr(A † A)Tr(B † B) .

‡ The necessity of this change of perspective in dealing with systems of identical particles, with emphasis on observables rather than states, has already been remarked in[11]-[13].