PhD Guy Dabi Gab-Leyba
email: guy_dabi.gab-leyba@etu.udamail.fr

Bertrand Laporte
email: bertrand.laporte@udamail.fr

Dabi Guy

Gab-Leyba

Keywords: Government Take, Oil revenues, Oil tax, Monte Carlo simulation, Uncertainty, Chad. JEL Code : E27, E62, H21, N5

Concession Contracts (CC) and Production Sharing Contracts (PSC) have quite different implications for Government Take and the properties of the tax system, such as progressivity. In general, taxation via CC introduces significant distortions in activity, particularly due to the balance of royalties which tax production irrespective of the profitability of the project. So CC is normally regressive while PSC is normally progressive, because PSC taxation depends more directly on the profitability of the project. Chad has the distinction of having introduced PSC in the 2007 Chad oil code, while maintaining a royalty on production. Despite this feature, we show with a Cash Flow model and Monte Carlo simulations that the application of the 2007 oil code introduced more progressivity into taxation. This feature is particularly interesting in the current context of falling crude oil prices, because it maintains a favorable tax regime for exploration and exploitation by multinational oil companies. As a result, the Chad government should reactivate a counter-cyclical policy of oil revenue reserves when the crude oil price increases again.

Introduction

Industries such as those in the oil extraction sector can generate considerable income. Applying a tax system that captures a big part of the rent while providing incentives for exploration and development of new oil fields by International Oil Companies (IOCs) is therefore the priority of oil producer Developing Countries. There is a specific tax system in the oil sector because of the sector's special features, including massive investment and a necessarily complex technology controlled by IOCs [START_REF] Petermann | Mining and corruption[END_REF]. The issue of rent sharing, and therefore the tax system, is a major challenge for Chad, whose economy is heavily dependent on oil. In 2013 oil accounted for over 70% of Chad's tax revenues, 90% of total exports, and 30% of nominal GDP (IMF,Report No. 13/284).

The discovery of oil in Chad goes back more than 40 years, but political tensions and major logistical barriers to the export of oil delayed the exploitation of crude until 2003. From 2000, the World Bank support accelerated the construction of a pipeline. At a cost of USD 4.2 billion, this pipeline then represented the biggest private investment in Sub-Saharan Africa. Chad had proven reserves of crude of around 1.5 billion barrels at end 2013, which puts it in 10th place in Africa. After reaching a peak of 173,000 barrels per day in 2005, production steadily declined to 91,000 barrels per day in 2013.The inauguration of new fields in 2014 stopped the downward trend in oil production. These fields should stabilize the production level at around 88 million barrels in 2017 (IMF,Report No. 14/282). In the absence of new discoveries, production will decrease gradually and be completed by 2034.

A few years after the beginning of oil extraction, a major disagreement occurred between the government and the IOC consortium1 about the oil revenue-sharing agreements. The 19882 Agreement covers the extraction of oil from three of the original Doba fields3 . The Petroleum Code of 1962, which gives the outlines of the contractual relationship between Chad and the IOCs, is based on CC.

The two main tax instruments are royalties and corporate tax.

The low Government Take4 is due to the taxation arrangements which were designed to attract IOCs because of the low quality of Chad crude oil, high transport costs (no coastline, construction of a pipeline), and generally high operating risks [START_REF] Blaise | Fiscalité pétrolière au sud du Sahara : la répartition des rentes[END_REF]. The government's inexperience in negotiations with the IOCs is another reason which explains that this share was deemed unfair by the government (Margonelli, 2007).

In 2004 the government tried to renegotiate the tax system of the 1998 Agreement. Despite the 1.6 million USD in legal aid provided by the World Bank, the renegotiation has not resolved the disagreement on Government Take, because of a non-retroactivity clause which locks the key provisions of the 1988 Agreement [START_REF] Gary | Le pétrole tchadien : miracle ou mirage? Suivre l'argent au dernierné des pétro-États d'Afrique[END_REF].This disagreement crystallized following the rise in crude oil prices since 2004, which allowed a rapid return on investment5 . To modify the terms of sharing of oil revenues, the new 2007 oil code substitutes the Production Sharing Contract (PSC) for the CC.

There is no "standard" CC or PSC. These are different instruments and methods of taxation that impact differently in the way they approach oil revenue sharing. Both types of contract can therefore, in principle, enable the Government to achieve its objectives of oil revenue sharing. Given the importance of oil to the Chad economy and the need for further exploration activities and field development, despite the unfavorable economic climate and the fall in oil prices, this paper compares the two kind of contracts used in Chad mainly through Government Take and the incentive nature of the tax design. The rest of the paper is organized in three sections. The first section presents the oil tax design and its evolution. The second part presents the literature and describes the methodology and data used in our model. The third gives the results and analyzes the changes between the two agreements.

The tax regime design for the upstream oil sector in

Concession Contract Agreements

For the operations of the Doba basin, three legal documents bind the IOC Consortium, Chad, and Cameroon. Firstly the establishment of the two agreements for the operating companies for; a. the Cameroon side Chad pipeline downstream, and b. "research convention, exploitation and transportation of hydrocarbons" upstream, both dating from 1988 and which were subject to later revision (amendment No. 2, March 1997), these two documents set out the modalities for the management of the pipeline (royalties, taxation and law, litigation between the parties, etc.). The third between Chad and the Consortium only, establishes the taxes levied on operations.

The Convention of 1988, supplemented by the 2004 agreement, based on the 1962 oil code enshrines the concession system as the sharing mode for oil revenues. The tax provisions of the Code concern mainly the application of a royalty and a net profit tax (corporate tax). The royalty rate is fixed in the 1962 code at 12.5%. It was increased in 2004 to 14.25% at the signing of the agreement governing the new fields of Timbré and Maikeri. The net profit tax rate depends on the field. Apart from the Nya and Moundouli fields to which is applied a single rate of 50%, the other fields are subject to a rate that varies with the price of crude and an R factor, calculated as the ratio between revenue accumulated net after-tax upstream of the investor and the investors' cumulative upstream investment costs. These agreements allow for several exemptions -the contractual minimum tax (IMF), patent fees, land tax, Value Added Tax (VAT), duties and taxes on petroleum products, vehicle taxes, levies on dividend distributions, and technical assistance. Lastly, the 1988 Convention contains a stability clause (Article 34) which "fixes" the terms of sharing of oil revenues.

Production Sharing Contract Agreements

The 2007 oil code and the application laws that followed were initiated to increase the share of oil revenues collected by the Government. The 2007 oil code maintains a royalty on production at rates ranging from 14.25 to 16.5% and introduces a public capital stake of the Societé des Hydrocarbures du Tchad (SHT) 6 in the capital of the private operator, which may not exceed 20%. The maximum public capital stake was set in 2014 at 25% in future farm projects following the renegotiation of exploration contracts with the China National Petroleum Corporation International (CNPCIC). This state participation is made in return for payment without interest in proportion to the costs of the search and development operations of oil fields (Paragraphs 4a and 4b of the Petroleum Code's Article 14). 6 A public company representing the interests of the Government.

The production sharing contract between SHT and the private operator shares net crude oil production that is the crude oil production less the amount of crude oil which covers production costs (cost oil).

The code states that these costs cannot exceed 70% of crude oil production minus the royalties levied on all contractual areas. The remainder of the crude oil production, the profit oil, is split between the private operator and the SHT according to an R factor (Table 2). The numerator is the amount of sales accumulated since the beginning of the operation until the last quarter before the quarter for the sale, less accumulated operating costs (operational costs and investments) in the same period. The denominator is the sum of the accumulated costs of development operations and costs of exploratory activities conducted. A statistical royalty is also levied on the value of oil exports. Its tax base is the value of exported crude oil, which is slightly different from the production of crude oil due to the existence of an amount taken from the pipeline to operate the pumping centers. The statistical royalty rate is 1% for the 1998 agreement, and 1.5% for the 2004 agreement. The 2007 oil code introduced two types of bonuses. The Signature Bonuses whose amount is not fixed, and Attribution Bonuses for an exclusive operating license whose amount is set at $ 2 million. These two types of bonuses are considered as an unrecoverable oil cost.

Methodology and data

The evaluation of tax design for natural resources has been the subject of an abundant literature both theoretically and empirically. [START_REF] Smith | Modeling the Impact of Petroleum Exploration and Development[END_REF] provides a broad review of work highlighting the wide variety of methods and models. [START_REF] Andon | Comparing petroleum fiscal regimes under oil uncertainty[END_REF]Roberts (2006), following Lund (1992), use the approach known as "contingent claim analysis" 9 to compare 5 tax systems 10 . In a context of uncertainty about prices and using Monte Carlo simulation techniques, they classify these tax regimes according to the economic distortions they cause. [START_REF] Diderik | Petroleum taxation under uncertainty: contingent claims analysis with an application to Norway[END_REF], in an application to the Norway tax regime shows that the application of concession contracts lead to significant distortions for the Norwegian economy.

Many authors use more empirical methods based on the development of scenarios. A basic scenario defines the volume of reserves, the scope and the investment period, the number of wells, the drilling rate success, the intensity of development, the initial production rate (or decline), the variable operating costs, etc. Based on assumptions about oil prices and the main taxes, it is then possible to calculate the cash flow and the Government Take over the project lifetime, for different scenarios.

The Peaceman model ( 1977) is based on the simulation of the crude extraction rate of gas and water from a reservoir. Reservoir characteristics are obtained by exploiting the geophysical modeling and three-dimensional geology to capture heterogeneous reservoir features. By simulating the physical and financial consequences of an alternative drilling program, this approach allows for optimization of the value of the resource. The investor's behavior can be studied after the introduction of a tax which would change the structure of net cash flows, regardless of the anticipated development plan of the field. [START_REF] Uhler | The Rate of Petroleum Exploration and Extraction[END_REF] built an integrated model of the exploration and extraction process which, when applied in its simplest form to a homogeneous reservoir leads directly to a model of exponential decline in which the rate of production of each field decreases in a constant rate each period. This approach also takes into account the pressure injection operations that can slow the rate of decline and increase the volume of recoverable reserves. [START_REF] Jacoby | The Effects of Taxes and Price Regulation on Offshore Gas[END_REF], Helmioskoui and others (1992), and [START_REF] Rao | An Integrated Modeling Framework for Exploration and Extraction of Petroleum Resources[END_REF] have refined the model by incorporating other variables such as the heterogeneity of reservoirs and pressure.

There is a wide variety of papers on how taxation affects different aspects of the development of oil projects. These papers include a variety of methods to analyze the behavior of IOCs. Like Poterba (2010), Smith (2012 11 ) considered the ideal approach would be to model a comprehensive tax system 9 This approach allows consideration of the cash flows as a nonlinear function of the value of the underlying asset. based on a robust production model of an extractive industry. [START_REF] Smith | Modeling the Impact of Petroleum Exploration and Development[END_REF] proposes a modeling technique which specifically permits the analysis of how IOCs would adjust the intensity of exploratory activities, the timing and intensity of the initial development of the fields, the timing and intensity of enhanced recovery, and the eventual abandonment of fields in unfavorable conditions due to lower prices. This is a simple exploration and development model for studying the performance of an alternative tax system, and to identify potential distortions. The performance of a natural resource tax system depends on (i) its ability to increase government revenues; (ii) potential distortions that reduce the value of the after tax cash flows; (iii) the allocation of risk between the Government and the IOC.

To compare the two different agreements in Chad, we use a cash-flow model applied to real data from the Doba project, covering its technical features and cost structure. The analysis focuses on the Government Take, the progressivity in taxation, and the Internal Rate of Return (IRR) for the private operators of the two oil agreements, after making a series of assumptions about changes in oil prices and the rate of extraction (Monte Carlo simulations).

Production Profile

During the initial extraction, the extraction rate is determined by natural forces in the reservoir (pressure, permeability of reservoir rocks, viscosity of crude trapped, etc.) [START_REF] Smith | Modeling the Impact of Petroleum Exploration and Development[END_REF]. Production (Qt) then decreases from the initial production level (Q 0 ) at a fixed rate (a) during that initial production phase:

𝑄 𝑡 = 𝑄 0 𝑒 -𝑎𝑡 𝑡 ≥ 0 (1)
The production falls due to the drop in pressure in the reservoir as the volume of crude decreases [START_REF] Uhler | The Rate of Petroleum Exploration and Extraction[END_REF]. The depletion rate, and therefore the rate of decline, is determined by the intensity of the operator's development effort (number of wells drilled in the reservoir, drilling platforms, and water injectors).

The volume of primary reserves (R) is given by the following formula:

𝑅 = ∫ 𝑄 𝑡 ∞ 0 𝑑𝑡 = 𝑄 0 𝑎 (2)
According to the equation ( 2), 𝑄 0 = 𝑎 × 𝑅 .The rate of decline (a) and the extraction rate are identical.

It follows from this equation that the volume of the remaining reserves in the reservoir (R t ) at time t is given by:

𝑅 𝑡 = 𝑄 𝑡 𝑎 = 𝑅𝑒 -𝑎𝑡 (3)
The first year of production is given by the formula: Suslick (2003); [START_REF] Johnston | International Exploration Economics, Risk, and Contract Analysis[END_REF]; [START_REF] Tordo | Fiscal Systems for Hydrocarbons: Design Issues[END_REF]; Johnston, Johnston and Rogers (2008);and Daniel et al (2010).

q 1 =a×R, (4) and production for the following years by the formula: q t+1 = (1-a)×q t.

(5)

The cash-flow model

The Net Present Value (NPV) of the project before government tax can be calculated from the following formula:

NPV = ∑ 𝑍 𝑡 -𝐶 𝑡 -𝐾 𝑡 (1+𝑖 * ) 𝑛 𝑇 𝑡=0 ( 6
)
𝑍 𝑡 is the expected turnover from the sale of crude oil, C t the unit operational cost expressed in $/b, K t the capital cost for the initial development of operations and the maintenance of production.

The internal rate of return (IRR) is used to assess the profitability of the project. The OIC compares it to the minimum threshold expected to develop the business (Macey 1998). In this sector, OICs develop the fields only if feasibility studies show an IRR greater than or equal to 100% at the latest during the seventh year of production. That is explained by the fact this big return on investment is needed to minimize project risks. It is therefore a decision tool for investment. The before tax IRR of a project is obtained from the following equation:

NPV = 0 = ∑ 𝑍 𝑡 -𝐶 𝑡 -𝐾 𝑡 (1+𝑖) 𝑛 𝑇 𝑡=0 (7)
It corresponds to the yield or compound interest rate that cancels the NPV.

To estimate the tax rent generated during all the life of a project, we used the Average Effective Tax Rate (AETR) which is the ratio of the sum of all the parts of the rent collected under the various provisions of taxation and the before NPV of the project. It is given by:

𝐴𝐸𝑇𝑅 = ∑ 𝐹 𝑡 𝑇 𝑡=0 ∑ 𝑍 𝑡 -𝐶 𝑡 𝑇 𝑡=0 -𝐾 𝑡 (8)
F t represents the sum of government tax levies.

Oil price projection and Monte Carlo simulation

The forecasting of future crude prices is an important step in the evaluation of oil tax regimes in a context of uncertainty. The high volatility of crude prices does not permit to detect a long-term relationship describing the evolution of the price.

That why in this paper, the Brent oil price which constitutes the reference for the Doba oil project is forecasted on the basis of a Geometric Brownian Model (GBM). Lund (1992), and Blake and Roberts (2006) used the same approach since its projections have the advantage of being more realistic than the Auto Regressive models (AR). Our model is derived from the model developed by [START_REF] Hull | Options, Futures and Other Derivative Securities[END_REF] whose functional form is:

𝑃(𝑇) = 𝑃(0)exp{[𝑟 -𝛼 -( µ 2
2 )] 𝑇 + µ𝛽√𝑇}(9)

Where T represents time at the end of each tax cycle, P (T) and P (0) the price of Brent at the period T and 0, 𝛼 the net marginal convenience yield from storage (convenience yield less marginal storage costs), and r the nominal risk free interest rate.

The standard deviation (µ) and the average annual price of Brent (β) over the 1977 to 2013 period are 29.9% and $ 39.75 per barrel respectively. The nominal rate of interest without risk is approximated as in Blake and Roberts ( 2006) by the rate of return of United States bonds for the 30 years corresponding to the life of the project. The data used to generate the price projections are shown in the table below. Given the extreme volatility of crude oil prices, including the sudden climb from 2005 when prices rose from $ 54 to more than $ 111 per barrel in 2012, our simulations are based on 5 price scenarios with 10,000 simulations per scenario. This Monte Carlo approach takes into account the uncertainty of the variables that influence the project profitability by using probability distributions. The PSC applied in the Doba Project includes the following tax terms : We ignored the signature bonuses, statistical royalty and surface fees, taxes on wages, and license fees for which revenues can be considered as negligible.



Results and analysis

The sharing of oil revenues can only be considered if the tax system has enough incentive for investors to engage in the search for oil fields and in their operation. The tax system of the various contracts is evaluated in three dimensions: the share of oil revenues collected by the government, the government revenue, and the incentive nature of the tax system through the rate of IRR for the investor and the progressivity in taxation. 12 It corresponds to the total development expenditure weighted by the size of the project field reserves. The same treatment is used for exploratory investments. 13 The unit operating cost is an average unit cost over the period 2003-2008. 14 This is an average over the period 2003-2008. It is relatively high according to the consortium, due to an increased presence of water in reservoirs requiring the drilling of new wells to maintain the production level.

The analysis of the results initially covers the basic scenario in which a Brent price of 50 USD per barrel is used, followed by a sensitivity analysis with three other price scenarios. Based on 10,000 simulations, the average NPV of the project is USD 5.9 billion, and the average IRR before tax is 105.5% for the base case.

Analysis of tax regimes under the base case scenario

The Average Effective Tax Rate (AETR), which represents the share of oil revenues collected by the government, is relatively higher for the CC (74.5%) than for the PSC (70.5%). The government revenues are lower under the PSC than the CC. The After-tax IRR is lower under the CC (24.1%) than the PSC (40.4%) (Table 5).

Analysis of other price scenarios

By taking into consideration various price scenarios, we can assess the progressivity of the two tax systems. The AETR of the CC decreases as the price rises. The AETR of the PSC rises as the price rises. When the oil price increases from $ 50 to $ 110, the AETR decreases from 74.5% to 66.9% with CC, whereas it rises from 70.5 to 74.7% with the PSC. The state's revenues are greater for the PSC than the CC, except for low oil prices (Table 6). The gap between the before, and after IRR increases sharply when the Brent price increases. In absolute terms, the government revenue increases sharply with the increase in the price of Brent, but in a smaller proportion with the PSC than with the CC.

Chart 4: Evolution of before tax IRR and after tax IRR

Conclusions

There is no real "standard" CC and PSC. Each country defines the tools and methods of taxation it intends to apply through these contracts. The instruments and modalities of the contracts determine the properties of the tax system. In general, taxation instruments accompanying concession contracts introduce significant distortions in activity, particularly due to the effect of royalties as a taxation instrument, because royalties tax production activity irrespective of the profitability of the project. So CC taxation is often regressive. The "pure" PSC is normally progressive because oil revenues are mainly based on profits tax, often with a rate which varies depending on the profitability of the project.

Chad has the specificity of having introduced into its legislation the PSC while maintaining royalties.

Although deductible on profit, this royalty is contrary to the expected progressivity of a tax system based on a PSC. Despite this feature, the application of Chad oil code 2007 allows the government to capture a larger share of oil revenues, compared to the 1988 agreement, but also a gradually increasing share as the profitability of the oil field improves. This feature is particularly interesting in the current context of falling crude oil prices, because it balances risks between State and Consortium, and maintains a favorable tax regime for exploration and exploitation by IOCs. As a result, the Chadian government should reactivate a counter-cyclical policy of oil revenue reserves when the crude oil price increases again.

.

  Chad: Concession Contract versus Production Sharing Contract Chad's oil tax regime is defined by its 2007 oil code (Act No. 006/PR/2007 of 2 February 2007), supplemented by Decree No. 796/PR/PM/MPE/2010 which lays down detailed rules, and Ordinance No. 001/PR/2010 approving the PSC. If these texts constitute the institutional framework for reference in legal and tax matters for the sector, the CC signed before 2007, under the Ordinance No.007/PC/TP/MH February 3 rd 1962, remains in effect. These are primarily the agreements signed in December 1998 and May 2004 concerning the operation of the Doba field, the agreement signed in February 1999 confirmed in 2012 by the CNPCIC, and the agreement signed in January 2006 confirmed in 2012 by OPIC Africa (Taiwanese national oil company).

10

  These regimes include those of the Province of Alberta in Canada which consists of a Corporate Income Tax and a Royalty; the pre-2003 regime of Papua New Guinea (PNG) based on the traditional Rate of Return (RR); the joint project of Sao Tome and Principe & Nigeria based on a PSC (SNDJZ); the Tanzanian hybrid system based on RR /PSC; and the PSC regime of Trinidad and Tobago.

Chart 1 :

 1 Changes in Brent oil price ($/barrel) Source: British Petroleum. Box 1: Basic parameters of the oil fields and tax systems used in the model  Proven reserves (Millions (Millions of Barrels) : 588  Rate of extraction (%) : 7.2  Development investment level of the field (USD MM) : 942.42 12  Exploratory cost ($ MM) : 435.5  Unit Operating Cost per barrel, OPEX ($/b) : 5 13  Unit Capital Cost, CAPEX ($/b) : 11 14  Discount rate r (%) : 8  Base case Brent price per barrel ($) : 50  Discount of Doba crude oil (%) : 10  Transportation Cost per barrel ($/b) : 10 Two kinds of tax regimes are analyzed in the model :

Table 1 :

 1 Rate of Petroleum net profit tax

	Crude Oil Reference Price US$1995/barrel	R ratio<1.75	1.75<R≤2.50	R≥2.50
	≤18$	40%	50%	65%
	>18$ et ≤ 19$	45%	52,5%	65%
	>19$ et ≤ 20$	50%	55%	65%
	>20$ et ≤ 21$	55%	57,5%	65%
	>21$ et ≤ 22$	55%	57,5%	65%
	> 22.00$	60%	60%	65%
	Source: Convention of 1988 for three originals fields (Komé, Miandoum and Bolobo)

Table 2

 2

	: Profit Oil distribution scheme
	R Factor	<2.25	2.25<R<3	>3
	IOC's share in the Profit Oil	60%	50%	40%
	OPC's share in the Profit Oil	40%	50%	60%
	Source: PSC between Chad and PETROCHAD MANGARA.

Table 3 :

 3 List of oil exploitation in Chad and nature of contract

	companies in production
		Operator	Bloc	Partner		Nature of contract	Date of signature	International Oil
								Companies
	1	Esso E&P CHAD (40%)	Chari Bolobo	Petronas Caligari	(35%) and	Concession Contract	December 1988 and June	EXXONMOBIL
			Chari Komé	Chevron Petroleum Chad (25%)		2004	(United States)
			Chari Maikeri
			Chari Timbré
			Chari Moundouli
			Chari Miandoum
			Chari Nya
	2	CNPCI (100%)	H (Rônier)			Concession Contract	February 1999	CNPC (China)
			H (Mimosa)
			H (Madiago)
	3	GLENCORE ENERGY (75%)	DOB	SHT (25%)		Production Sharing Contract		Glencore-Xstrata
			DOI					(Switzerland)
	Companies in exploration
	1	CNPCIC (100%)	H (Bongor Est)			Concession Contract	February 1999	CNPC (China)
			H (BongorOuest)
			H (Chari Ouest)
			H (Lac Tchad)
			H (Madiago Est)
			H (MadiagoOuest)
			H (Sédigui)
	2	OPIC (70%)	Chari Ouest III (50%)	SHT (30%)		Concession Contract	January 2006	CPC (Taiwan)
			Chari Sud II
			Lac Tchad I
	3	ERHC ENERGY (100%)	BDS-2008			Production Sharing Contract		ERHC Energy Inc.
		ERHC ENERGY (100%)	Chari Ouest III (50%)					(United States)
		ERHC ENERGY (100%)	Manga
	4	GLOBAL PETROLEUM (100%)	Djado III			Production Sharing Contract		Nigeria
			DOE
			DOF
			DOG
			Largeau V
	5	GLENCORE ENERGY (100%)	Borogop			Production Sharing Contract		Glencore-Xstrata
		GLENCORE ENERGY(100%)	Chari Est Doséo					(Switzerland)
		GLENCORE ENERGY(100%)	DOH
	7	PETRA BV (100%)	Erdis			Production Sharing Contract		PETRA ENERGIA,

7 Contract cancelled in March (ETIE-Chad, 2012 Report, Fair Links.) 8 Contract cancelled in October 2011

(ETIE-Chad, 2012 Report, Fair Links.)

Table 4 :

 4 Parameters used in GBM modelParameters

		Values
	µ	29.9%
	β	$ 39.75
	P(0)	50
	α	4%
	r	4.75%

  Production Royalty (RP) with a rate of 12.5%  Corporate Tax (CT) with a rate between 40 and 65% and according to the Brent oil price.The Production Sharing Contract is based on the agreement concluded in 2011 for the exploitation of Badilla and Mangara fields located in the Dosseo basin. It involves the following tax terms :  Cost Oil must not exceed 70%  Profit Oil is shared between the OIC and the government through the national oil company SHT.The share is done according to a share rate depending on the R factor. The sharing of Profit Oil between the two partners is done as follows: o If R< 2.25 then the government share of Profit Oil is 40% o If 2.25<R< then the government share of Profit Oil is 50% o IF R>3 then the government share of the Profit Oil is 60%  Public participation in the capital of the consortium is 25%

Table 5

 5

	: Economic results in the base case
		AETR	IOC share of	After-tax IRR	Government	IOC revenues
		(%)	oil revenues	(%)	revenues ($	(($ MM)
			(%)		MM)
	Concession Contract	74.5%	25.5%	24.1%	4,605.4	1,304.4
	Production Sharing Contract	70.%	29.5%	40.4%	4,186.7	1,723.1

Table 6 :

 6 Economic performance under the other price scenarios Price ScenarioOil price :

$/b 110 Oil price : $/b 90 Oil price : $/b 70 Oil price : $/b 50

	Chart 3: Evolution of the after tax project profitability
	Before tax IRR (%)	525%	358%	220%	105%
	NPV (Million USD)	22,868.0	17,215.3	11,547.6	5,909.8
	AETR CC (%)	66.9%	67.7%	69.3%	74.5%
	Share of the IOCCC (%)	33%	32%	31%	25%
	After tax IRR (%)	159%	107%	63%	24%
	Government revenues CC (Million USD)	15,497.5	11,860.6	8,215.8	4,605.4
	IOC revenues CC (Million USD)	7,370.4	5,354.7	3,331.8	1,304.4
	AETRCPP (%)	74.7%	74.5%	73.9%	70.5%
	IOC share of the PSC (%)	25%	25%	26%	29%
	Before tax IRR (%) PSC (%)	168%	129%	88%	40%
	Government revenues PSC (Million USD)	17,032.0	12,769.3	8,451.3	4,186.7
	IOC revenues PSC (Million USD)	5,836.0	4,445.9	3,096.3	1,723.1

i Appendix 1: Mathematical formulae of the tax provisions used in the model

  BTCF 𝑡 = 𝑃𝐵 𝑡 𝑄 𝑡 -𝐶𝐴𝑃𝐸𝑋 𝑡 -𝑂𝑃𝐸𝑋 𝑡 -𝑇𝐴𝑋 𝑡 Where 2. TAX 𝑡 = 𝑅𝑃 𝑡 + 𝐼𝑆 𝑡 3. PR 𝑡 = 𝑟 1 (𝑃𝐵 𝑡 𝑄 𝑡 ) 4. CT 𝑡 = 𝑟 2 (𝑃𝐵 𝑡 𝑄 𝑡 -PR 𝑡 -𝐶𝐴𝑃𝐸𝑋 𝑡 -𝑂𝑃𝐸𝑋 𝑡 -𝐷𝐴𝑀 𝑡 ) Production Sharing Contract 5. PR 𝑡 = 𝑟 3 (𝑃𝐵 𝑡 𝑄 𝑡 ) 6. CS 𝑡 = 70%(𝑃𝐵 𝑡 𝑄 𝑡 -PR 𝑡 ) 7. CO 𝑡 = 𝐶𝐴𝑃𝐸𝑋 𝑡 + 𝑂𝑃𝐸𝑋 𝑡 𝑖𝑓 CO 𝑡 < CS 𝑡 ; Otherwise 70%(𝑃𝐵 𝑡 𝑄 𝑡 -PR 𝑡 ) for the first period 8. CO 𝑡 = 𝐶𝐴𝑃𝐸𝑋 𝑡 + 𝑂𝑃𝐸𝑋 𝑡 + (CO 𝑡-1 -CS 𝑡-1 ) 𝑖𝑓 CO 𝑡 < CS 𝑡 ; Otherwise 70% (𝑃𝐵 𝑡 𝑄 𝑡 -PR 𝑡 ) for the other periods 9. R = ∑ 𝐶𝐴𝑃𝐸𝑋 𝑡 + 𝑂𝑃𝐸𝑋 𝑡 + 𝑃𝐵 𝑡 𝑄 𝑡

	OPEX t =oQ t
	CAPEX t =cQ t
	PB t =(1-D)P t
	Concession Contract
	1. 𝑛 𝑡=0
			∑ 𝐶𝐴𝑃𝐸𝑋 𝑡 𝑛 𝑡=0
		𝑅 < 2.25	0 ⇒ 𝑟 4 = 40%
	10. If (	2.25 < 𝑅 < 3	0 ⇒ 𝑟 4 = 50%	)
		𝑅 > 3	0 ⇒ 𝑟 4 = 60%

0 11. PO 𝑡 = 𝐶𝐴 𝑡 -RP 𝑡 -CO 𝑡

Initially Chevron, Exxon Mobil, and Petronas.

The first version of this agreement of the CC kind was signed on December 19 th 1988. This version was amended on May 19 th 1993, March 12 th 1997 and finally October 16 th 1999.

Komé, Miandoum and Bolobo.

The Government Take for a given field, a concession or a basin is the ratio of the sum of the undiscounted net cash flow accruing to the State over the entire commercial life of the deposit based on the total net Cash flow generated by oil operations (including related research investment deposit, concession or basin).(Leenhardt B., 2005).

According to the 2012 Exxon Mobil Annual Report, the rise in crude prices starting in 2004 allowed a rapid return on investments, and early payment of Corporate Tax in 2006 (Project Update No. 33).

Coefficient of variation. It corresponds to the ratio of the standard deviation to the mean

This work was supported by the LABEX IDGM+ (ANR-10-LABX-14-01) within the program "Investissements d'Avenir" operated by the French National Research Agency (ANR). http://www.cerdi.org/

Data and sources

The tax system based on the CC is regressive while that based on the PSC is progressive (Chart 3).

The AETR of the CC curve clearly decreases whereas that for CPP increases slightly before stabilizing.

Chart 2 : Evolution of the progressivity of the two tax regimes

The After-tax IRR of the PSC is always higher than that of the CC, although the two curves converge at high oil prices (Chart 4). The two contracts are equivalent for the investor at high prices but diverge significantly at low prices. 100 simulations at each dollar level of price between $ 57 and $ 61 per barrel were conducted to determine the price at which government revenues and the AETR of the PSC become higher than those of the CC. For this price level, government revenues are around 5.8 billion of $ for an after-tax IRR of 36.4%.


C Ducoin

Ph Chomaz

F Gulminelli

Isospin fractionation : equilibrium versus spinodal decomposition

Keywords: numbers: 24.10.-i, 24.10.Pa, 24.60.Ky, 25.70.Mn, 64.70.Fx, 65.40.Gr

This paper focuses on the isospin properties of the asymmetric nuclear-matter liquid-gas phase transition analyzed in a mean-field approach, using Skyrme effective interactions. We compare two different mechanisms of phase separation for low-density matter: equilibrium and spinodal decomposition. The isospin properties of the phases are deduced from the free-energy curvature, which contains information both on the average isospin content and on the system fluctuations. Some implications on experimentally accessible isospin observables are presented.

I. INTRODUCTION

The nuclear interaction is a mean-range attractive, short-range repulsive force. Because of this, nuclear matter is expected to present a phase transition of the liquid-gas type [1,2].

The nuclear-matter phase diagram has been studied in detail in the case of symmetric matter, i.e. a system with the same number of protons and neutrons [3]. The case of asymmetric matter involves an additional degree of freedom, isospin. As a result, its study is more complex and allows to observe specific effects such as the phenomenon of isospin fractionation [4]. The knowledge of these properties is of great interest for the physics of neutron-rich systems. Such systems have become the subject of many investigations in the recent nuclearphysics literature. Exotic-nucleus transport properties and multifragmentation are expected to provide constraints on the nuclear equation of state, in particular on the density dependence of the symmetry energy [5,6,7]. In astrophysics, both neutron-star structure and supernova dynamics are influenced by phase transitions of neutron-rich matter in a large interval of temperatures and densities [8,9]. All these different phenomena involve excited matter at baryon densities lower than normal nuclear matter density : this corresponds in the phase diagram to a region of instability with respect to phase separation. In the case of neutron-star crusts the system lifetime is such that one can safely modelize phase separation in terms of equilibrium. However this is not the case for exploding supernova cores, and even less for heavy-ion collisions. For such systems with a finite lifetime, the out-of-equilibrium phenomenon of spinodal decomposition may be the leading mechanism for phase separation [10].

For the work presented in this paper, we have studied nuclear matter in a mean-field approach, using Skyrme effective interactions. Three different parameterizations are used.

The more recent one, Sly230a [11], is used as a reference since it has been fitted to exotic nuclei and neutron-matter properties, which makes it particularly well-adapted to the description of neutron-rich matter. Two others are used as a comparison : SGII [12], which has been shown to reproduce isospin effects studied through giant dipole resonances [13],

and SIII [14], which is one of the first Skyrme parametrizations and is not constrained for asymmetric matter. Contrary to Sly230a and SGII, SIII is certainly not a realistic parameterization; we have included it in our comparison since it is still used in some transport codes because of its simplicity [5].

In section II, we present the phase diagram of nuclear matter, which is a representation of phase separation at equilibrium. Then the phase properties are considered, comparing the case of equilibrium and spinodal decomposition. In section III, we study the direction of phase separation, leading to the phenomenon of isospin fractionation. The consequences of our findings on the possible extraction of the symmetry-energy coefficient from isoscaling data [6,7] are explored in section IV. In section V we discuss fluctuations in the system isotopic composition.

II. PHASE DIAGRAM

Gibbs statistical equilibrium is defined by the maximization of the system Shannon entropy S = -T r( Dlog D), where D is the many-body density matrix of the system [START_REF] Balian | From Microphysics to Macrophysics[END_REF]. Let us consider S represented for a homogeneous system in the space of observables. If S presents a convexity in this representation, it can be maximized by phase mixing, which corresponds to a linear interpolation between two points of the space of observables. The coexistence region is then defined as the ensemble of points for which entropy can be enhanced by such an interpolation. If the system is located at a given point inside coexistence region, there is a unique decomposition in a couple of two other points that gives the best maximization of entropy : they define the two phases at equilibrium. This couple is determined by Gibbs construction [START_REF] Huang | Statistical Mechanics[END_REF], which consists in building the concave envelope of the entropy surface by linear interpolation in the space of observables : phases at equilibrium then belong to the same tangential plane [9,[START_REF] Muller | [END_REF]. This geometrical condition results in the well-known equality of all intensive parameters.

Let us consider a system occupying a volume V described by a set of observable densities {a k } = { Âk /V } controlled by the associated intensive parameters {λ k }, related to the entropy density s = S/V by :

λ l = ∂ a l s({a k }) (1)
The pressure P is related to the entropy by :

βP = s - k λ k a k (2)
Let us now consider two points of the space of observables {a k } (1) and {a k } (2) . Entropy has the same tangent plane at these two points if {λ k } (1) = {λ k } (2) (equality of plane slopes) and P (1) = P (2) (equality of plane origins) : these are Gibbs equilibrium conditions on the intensive parameters.

In the case of asymmetric nuclear matter, the relevant observables are neutron, proton and energy densities, {a k } = {ρ n , ρ p , e}; the associated intensive parameters are {λ k } = {-βµ n , -βµ p , β} where µ q (q = n, p) are the chemical potentials and β is the inverse temperature. To obtain the phase diagram of nuclear matter, expressions linking observables and the associated intensive parameters (i.e. equations of states) are needed.

We have worked out such relations in the mean-field approach [18].

Let us first express the average energy density of homogeneous nuclear matter, e = Ĥ /V . With a Skyrme effective interaction, e is a functional of one-body densities only: particle densities ρ q and kinetic densities τ q = P 2 q / 2 V . We introduce isoscalar and isovector densities :

ρ = ρ n + ρ p , τ = τ n + τ p ρ 3 = ρ n -ρ p , τ 3 = τ n -τ p (3)
In the case of homogeneous, spin-saturated matter with no Coulomb interaction, the energy density can be written as :

e = 2 2m τ + C 0 ρ 2 + D 0 ρ 2 3 + C 3 ρ σ+2 + D 3 ρ σ ρ 2 3 + C ef f ρτ + D ef f ρ 3 τ 3 (4)
where coefficients C i an D i are linear combinations of the standard Skyrme parameters [11,18] C 0 = 3t 0 /8 D 0 = -t 0 (2x 0 + 1)/8

C 3 = t 3 /16 D 3 = -t 3 (2x 3 + 1)/48 C ef f = [3t 1 + t 2 (4x 2 + 5)]/ 16
D ef f = [t 2 (2x 2 + 1) -t 1 (2x 1 + 1)]/16 .
Mean-field approaches describe nucleons as independent particles. Each particle of type q is submitted to the mean field Ŵq defined by the relation δ Ĥ = T r( Ŵq δ ρq ), where δ ρq represents any variation of the one-body density matrix associated with the particle q, leading to the expression :

Ŵq = ∂e ∂τ q p2 2 + ∂e ∂ρ q = 1 2m * q p2 + U q (5)
The mean field is composed of a potential term U q = ∂ ρq e, and a kinetic term involving an effective mass m * q defined by 2 /2m * q = ∂ τq e. This determines individual energy levels ǫ i q = p 2 i 2m * q + U q with the usual quantification conditions on the momenta provided by the boundary conditions. In the grand-canonical approach, for a given temperature T = β -1 , these levels are occupied according to Fermi-Dirac statistics :

n q (p) = 1/ 1 + exp(β(p 2 /2m * q + U q -µ q ) (6)
where µ q is the chemical potential of the particles of type q .

All one-body densities are then given by Fermi integrals. In particular the independentparticle partition sum can be expressed in the same way as for an ideal gas:

ln Z 0 V = 2 q ∞ 0 ln(1 + e -β( p 2 2m * q -µ ′ q ) ) 4πp 2 h 3 dp = q 2 3m * q βτ q (7)
where µ ′ q = µ q -U q . At this point it is important to recall that the independent particle partition sum Z 0 should not be confused with the mean-field grand-canonical partition sum which is related to the mean-field pressure log Z = βP V and can be computed from the mean-field entropy using log

Z = S -β(< Ĥ > -µ p < Np > -µ n < Nn >).
For fixed values of ρ n and ρ p at a given temperature, chemical potentials are obtained by iteratively solving the self-consistent relations between ρ q , µ q and m * q . We can then compute the kinetic density and calculate the entropy as a function of densities. Indeed, the mean-field solution comes from a variational approximation to the entropy so the mean-field entropy S 0 provides the best independent-particle approximation to the exact entropy. The mean-field entropy can easily be derived from the independent-particle partition sum and the average of the mean-field potential Ŵ 0 so we get [19] :

S/V ≃ S 0 /V = ln Z 0 /V + β( Ŵ 0 /V -µ n ρ n -µ p ρ p ) (8)
which gives for a fixed β the entropy s = S/V as a function of neutron and proton densities.

As we have already mentioned, the construction of the concave envelope of s(e, ρ n , ρ p ) implies Gibbs conditions for each couple of phases at equilibrium :

β (1) = β (2) , µ (1)
n = µ (2) n , µ (1)
p = µ (2)
p and P (1) = P (2) . Let us now introduce the constrained entropy defined by

s c = -βf = s -βe. (9)
which corresponds to the opposite of the free energy f , divided by the temperature. Gibbs construction can then be performed by constructing the concave envelope of s c in (ρ n , ρ p ) plane at a fixed temperature T = β -1 . Indeed, β is imposed to be constant, and the definition of s c ensures that its tangent-plane slopes and altitude give respectively the chemical potentials and the pressure :

∂ ρq s c (ρ n , ρ p ) β = ∂ ρq s(e, ρ n , ρ p ) = -βµ q (10)
s c + β(µ n ρ n + µ p ρ p ) = s -β(e -µ n ρ n -µ p ρ p ) = βP (11)
The introduction of the constrained entropy permits then to reduce the dimensionality of the problem since we can work at constant temperature. The constrained entropy is represented in fig. 1, at a temperature of 10 MeV. In this figure, a slope has been subtracted in order to emphasize curvature properties: the surface represented is s ′ c = s c + βµ s ρ, where µ s is the chemical potential of symmetric matter at the transition. The constrained entropy is symmetric with respect to the axis ρ n = ρ p , reflecting the invariance of nuclear interaction with respect to neutron-proton exchange. We observe that the surface obtained presents a region of positive convexity : this is the case if we fix any temperature under the critical temperature of symmetric matter, which is, in the Sly230a parameterization, T c = 14.54MeV (for comparison, we get T c = 14.46MeV with SGII and 17.96MeV with SIII). The concave-envelope construction performed on this surface links couples of phases at equilibrium. The ensemble of these points forms the two branches of the coexistence curve, one at lower density (gas) and the other at higher density (liquid). They join at two critical points, where the difference between the two phases disappears.

The coexistence curve is the boundary of the coexistence region, inside which a system at equilibrium is decomposed into two phases situated on each branch. The region presenting a positive convexity is the spinodal region. In this region, infinitesimal density fluctuations may lead to an increase of entropy. The system is locally unstable.

III. DIRECTION OF PHASE SEPARATION

Let us consider a system whose proton fraction Z/A and total density ρ are such that it is situated inside instability region. This is illustrated on the right-hand side of fig. 2 for a neutron-rich system (Z/A = 0.3).

If equilibrium is reached, this system will undergo phase separation according to Gibbs construction. We can see that the corresponding phases in coexistence, represented as black dots on the coexistence border on the right part of fig. 2, do not belong to the line of constant proton fraction. We can remark that the liquid fraction is closer to symmetric matter than the gas fraction: this behavior is a consequence of the symmetry-energy minimization in the dense phase. This inequal repartition of isospin between the two phases is the well-known phenomenon of isospin fractionation [START_REF] Muller | [END_REF].

In order to quantify the phenomenon of isospin fractionation, we will study the associated directions in the density plane. The line of constant proton fraction defines a direction u Z/A .

Let us introduce the isoscalar direction u ρ , which corresponds to the axis of total density ρ = ρ n + ρ p . The equilibrium direction u eq is given by the straight line that joins the couple of phases in coexistence. Isospin fractionation is linked to the difference between u eq and u Z/A . The equilibrium direction being rotated towards u ρ , phase separation gives a liquid more asymmetric than the gas.

We now turn to study isospin fractionation if phase separation occurs out of equilibrium. Specifically, if the system is brought inside the spinodal region in a reaction too fast for global equilibrium to be achieved, its evolution will be driven by local instabilities instead of global Gibbs construction. This is spinodal decomposition [10] : the local properties of the constrained entropy curvature determine the development of a spinodal instability into phase separation. The information for this study is contained in the free-energy curvature matrix given by :

C =   ∂ 2 f /∂ρ 2 n ∂ 2 f /∂ρ n ∂ρ p ∂ 2 f /∂ρ p ∂ρ n ∂ 2 f /∂ρ 2 p   =   ∂µ n /∂ρ n ∂µ p /∂ρ n ∂µ n /∂ρ p ∂µ p /∂ρ p   ( 12
)
This matrix is defined for each point of the density plane. The lower eigen-value C < corresponds to the minimal curvature of f at the considered point. The spinodal border is the ensemble of points for which C < = 0. The instability region is defined by C < < 0 : in such points, there are directions of abnormal curvature, corresponding to a concave free energy. The eigen-vetctor u < associated with the lower eigen-value gives the direction of most negative curvature: this is the instability direction. It is represented by the double arrow on the example of fig. 2, where we can see that spinodal decomposition leads to a more pronounced fractionation than equilibrium, the denser phase getting closer to the symmetric matter. It is important to remark at this point that the instability direction does not give a direct quantitative measurement of the final fractionation after phase separation. For this to be true, the effect of the evolution of density fluctuations on the isospin composition has to be completely neglected. If this is certainly an approximation [27], it has been however

shown in explicit transport calculations [10] that the initial instability largely dominates the subsequent dynamics. The interest of employing this simplified approach as complementary to more sophisticated transport codes is that, within this approximation, the spinodal mechanism can be directly compared to phase separation at equilibrium. Indeed the equilibrium direction can also be expressed as the eigen-vector of a curvature matrix, like the instability direction. Let us define f eq as the free energy at equilibrium, i.e. after correction of the curvature anomaly by phase mixing according to Gibbs construction. This construction imposes a positive curvature at each point, and a zero curvature between two phases in coexistence. As a result, in the coexistence region the lower eigen-value of the f eq curvature matrix is zero, and the associated eigen-vector is what we have defined as the equilibrium direction, u eq . In the following, we will denote f 0 and f eq the free energy before and after Gibbs construction, their respective curvature eigen-modes being (C 0 < , u 0 < , C 0 > , u 0 > ) and (C eq < , u eq < , C eq > , u eq > ). A qualitative impression of how equilibrium and instability directions compare through the density plane can be gained from fig. 3 for the three Skyrme forces.

In order to follow the eigen-direction of the curvature matrix, we draw curves which are everywhere tangential to the curvature eigenvectors. Both directions of phase separation at equilibrium (Gibbs construction) and out of equilibrium (spinodal decomposition) remain quite close, and essentially isoscalar [5,20]. A comparable degree of isospin fractionation, with formation of a more symmetric liquid, is conserved all through the instability region.

The dominant feature is that the instability direction leads to more fractionation than equilibrium. In order to have a quantitative comparison of the different fractionations, it is convenient to express the directions of phase separation as slopes with respect to the isoscalar direction δρ 3 /δρ. This quantity gives a measure of the degree of fractionation associated with phase separation. It can in principle vary from zero (in the case of a purely isoscalar phase separation), to infinity (in the case of a purely isovector instability, corresponding to a separation of protons from neutrons into two phases at the same baryon density). In practice δρ 3 /δρ is always a small number [5,20]. This is easy to understand considering that i) proton and neutron kinetic energies favour equal-size Fermi spheres and ii) the nuclear interaction is more attractive in the proton-neutron than the proton-proton and neutron-neutron channels in any realistic nuclear physics model (T = 0 contribution), which favours the existence of a mixed isospin bound phase.

The evolution of δρ 3 /δρ (degree of fractionation) with ρ 3 /ρ (system asymmetry) is given by fig. 4 for three different total densities : below, at, and above the respective critical density ρ c (i.e the total density at critical points). In this representation, the direction of constant proton fraction is a diagonal straight line while the isoscalar direction is an horizontal line.

On this figure, equilibrium and instability direction are compared for the different Skyrme parameterizations.

Phase-separation evolutions verify some expected constraints. First, the isospin invariance of Skyrme interactions ensures that the graphs are symmetric with respect to ρ 3 = 0.

It also imposes that all directions become purely isoscalar in the case of symmetric matter, i.e. all curves must join at the origin (δρ 3 /δρ, ρ 3 /ρ) = (0, 0). Another constraint is given by the ending points of the curves. At these borders (|ρ 3 | = ρ), the vanishing of one species means that its chemical potential at finite temperature goes to -∞, implying an infinite curvature of the free energy in the corresponding direction. The minimal curvature u 0 < is then oriented along the axis of the remaining species, which is also the direction of constant

Z/A.
The first observation is that spinodal and phase-equilibrium directions are very close and lie between the isoscalar direction and the constant-proton-fraction one. This means that they lead to a rather similar fractionation. In order to get a deeper insight, let us consider the various densities in more details.

Let us consider the case of the axis at ρ = ρ c , for which the instability region ends at the critical points. Between these extremities, the instability direction induces more fractionation than equilibrium since it is closer to isoscalar direction. At the critical points, where the spinodal and coexistence curves are tangent, both directions of phase separation join to become tangential to the coexistence and spinodal borders. This is an additional constraint on phase-separation directions which, with the already discussed constraints, explains the similarities between the phase-separation directions at and out of equilibrium. This additional constraint does not apply if the total density we consider is not ρ c . We observe that for ρ < ρ c , a small inversion occurs around the extremities of the instability region. For ρ > ρ c instead, a phase separation out of equilibrium induces more fractionation for all isospin asymmetries. These observations stand for the three Skyrme forces that we have considered. The quantitative differences, which essentially distinguish SIII from the two more recent parameterizations, do not affect the qualitative features that have been described. This is partly due to the above discussed constraints.

IV. FRACTIONATION AND ISOCALING OBSERVABLES

In a variety of experimental situations, it has been observed that the isotopic-yield ratio for two reactions (denoted (1) and ( 2)) of different asymmetry approximatly follows an exponential form [6]:

R (1),(2) iso (N, Z) = Y (2) (N, Z)/Y (1) (N, Z) ∝ exp(α n N + α p Z). ( 13
)
This phenomenon is known under the name of isoscaling. In an equilibrium interpretation in the grand-canonical ensemble, the isoscaling parameters α q should be linked to the chemical potentials and temperature of the systems, α q = (µ

q -µ (2)
q )/T . It has been proposed to extract from the measured isoscaling parameters the symmetry-energy coefficient C sym through the approximate formula [7,23]:

α n = 4C sym ([Z/ A ] 2 (1) -[Z/ A ] 2 (2) )/T. ( 14
)
This relation, derived for the grand-canonical ensemble in the saddle-point approximation [7], holds in the hypothesis of equilibrium at the time of fragment formation. This means that to extract C sym (Z) from the measured isoscaling parameters, modifications of α n due to secondary decay have to be accounted for and the proton fraction of the most probable isobar for each Z should be known at the time of fragment formation .

Because of the small variation of α n with the atomic number, an average estimation of C sym can be obtained using for [Z/ A ] (i) the average proton content of the liquid fraction [7].

Since even this quantity is not experimentally known, a similar formula can be found in the literature, where, however, the proton fraction of the fragmenting source [Z 0 /A 0 ] (i) is used instead of [Z/ A ] (i) : this leads to a quantity we call C 0 sym [21,22,24]. This approximation is supported by the fact that in the zero-temperature limit the proton fraction of the fragments converges to the proton fraction of the source [22].

However, in general this does not need to be the case and the approximated C 0 sym may differ from the value C sym . In particular, both in the equilibrium hypothesis at finite tem-perature and in the case of dynamical fragment formation through spinodal decomposition, as we have seen, we can expect the phenomenon of fractionation: fragments should not have the same isospin ratio as the fragmenting source. This should affect the value of C sym obtained from the isoscaling coefficients.

Our equilibrium nuclear-matter approach does not allow to give quantitative predictions on the actual isospin content of fragments obtained out of the different phase-separation mechanisms; we can however make some qualitative order of magnitude estimations of the possible deformation that isospin fractionation may induce on symmetry-energy measurements.

Dynamical transport models [5,7] indicate that clusterization takes place at a relatively well-defined density inside the spinodal region determined by the collision dynamics and corresponding to the first development of instabilities, for which we take a typical value around ρ 0 /3. Spinodal decomposition is then a phenomenon fast enough for fragment formation dynamics to be dominated by the amplification of the most unstable modes at this density [10]. This density then determines the phase-separation direction.

The left part of fig. 5 displays, similar to fig. 4 above, the directions of phase equilibrium and spinodal decomposition as a function of Z/A at a typical temperature of T = 4MeV .

The case of two systems with global proton fractions corresponding to 112 Sn and 124 Sn is shown by the vertical lines. We can see that fractionation at equilibrium is negligible in agreement with ref. [22], and fragmentation of more exotic nuclei should be studied to explore this effect. However even for these stable systems, spinodal decomposition gives an important fractionation, and the effect would be amplified using more asymmetric systems.

Once the direction of phase separation is known, the final composition of the fragments depends only on their density at freeze-out, when they become separated. We call this density ρ sep . The average Z/A of the liquid fraction is represented in the central part of fig. 5 as a function of the density at separation. This density must be between the initial density and the saturation value, but its actual value is not known from fundamental principles. It is generally assumed that fragments at separation are at normal density ρ sep ≈ ρ 0 . Recent AMD calculations however seem to indicate that lower values ρ ≈ 0.1f m -3 may be more realistic [23,25].

The asymmetry coefficient C sym deduced through eq.( 14) is compared to its approximate value C 0 sym obtained with the proton fraction of the fragmenting source in the right part of fig. 5. Little change (few % at maximum) is observed in the case of equilibrium, again consistently with the predictions of the statistical model [22]. However, the stronger fractionation effect in the case of spinodal decomposition sensitively decreases, by up to 50% (or more), the approximate asymmetry coefficient. This implies that the low values reported for C sym in ref. [21] may be alternatively interpreted as a signal of isospin fractionation in a spinodal decomposition.

V. ISOSPIN FLUCTUATIONS

In section III, we have studied the direction of phase separation, corresponding to the direction of minimal curvature of the free energy u < . This direction can be qualitatively linked to the average isospin content of fragments issuing from the phase-separation mechanism.

As we have seen in section IV, some quantitative predictions can also be obtained if some extra hypotheses are made, namely that fragment-formation dynamics is dominated by an amplification of instabilities from well-defined initial density and temperature. Within the same line of reasoning, the other eigen-vector u > giving the direction of maximal curvature can be linked to fluctuations in the system composition, that we now turn to analyze.

A physical system characterized by a given value of the isoscalar density ρ at equilibrium presents a distribution of ρ 3 given by:

P β,µ,µ 3 (ρ 3 ) ∝ e V (s β,µ c (ρ 3 )+β µ 3 2 ρ 3 ) ( 15
)
where µ 3 = µ n -µ p is the isovector chemical potential, and µ = µ n + µ p the isoscalar chemical potential. The constrained entropy s β,µ c (ρ 3 ) is defined by the Legendre transform (cf eq.( 9)):

s β,µ c (ρ 3 ) = s(e, ρ, ρ 3 ) -βe + β µ 2 ρ (16)
In the case of an infinite system, the distribution ( 15) is characterized by a unique value,

P (ρ 3 ) ∝ δ(ρ 3 -ρ 0 3
), given by the maximum ρ 0 3 of the function s β,µ c (ρ 3 ) + β µ 3 2 ρ 3 . Indeed the correspondence between densities and chemical potentials is one-to-one for infinite systems since fluctuations go to zero at the thermodynamical limit. If conversely we consider systems of finite volume V, this distribution can be approximated by a gaussian if the system is big enough:

P β,µ,µ 3 (ρ 3 ) ∝ e - (ρ 3 -ρ 0 3 ) 2 2σ 2 ρ 3 (17)
where the distribution variance σ 2 ρ 3 =< ρ 2 3 > -< ρ 3 > 2 , which measures the isospin fluctuation, is directly linked to the entropy curvature. Indeed if the entropy has a finite curvature in ρ 3 direction, a development of s β,µ c (ρ 3 ) around ρ 0 3 gives:

s β,µ c (ρ 3 ) + β µ 3 2 ρ 3 = s β,µ c (ρ 0 3 ) + β µ 3 2 ρ 0 3 + (ρ 3 -ρ 0 3 ) 2 2 ∂ 2 s β,µ c ∂ρ 2 3 | ρ 0 3 (18)
so that

σ 2 ρ 3 = - 1 V ∂ 2 s β,µ c ∂ρ 2 3 -1 ρ 0 3 = 1 βV ∂ 2 f ∂ρ 2 3 -1 ρ 0 3 (19)
This is related to the symmetry-energy coefficient C sym defined from the free-energy curvature by ∂ 2 f /∂ρ 2 3 = 2 ρ C sym (ρ), which is independent of ρ 3 in the parabolic approximation [29] :

σ 2 ρ 3 = ρ 2βV C sym (ρ) (20)
The neglect of higher order terms in eq.( 18) is not fully justified for a genuinely finite system, where we expect deviations from the gaussian ansatz [START_REF] Muller | [END_REF]; it is however consistent with our mean-field approximation. Equations ( 17), (20) show that isospin distributions are directly linked to the symmetry properties of the effective interaction. This is why isotopic distributions measured in fragmentation reactions have been studied by many authors and their widths tentatively connected to the symmetry-energy coefficient C sym at subsaturation density [22,26].

Let us now consider a system at equilibrium inside the spinodal region at a point (ρ n , ρ p ).

Since this situation is unstable for an homogeneous system, it will separate into a liquid and a gas fraction. As we have discussed in section III, this separation follows the direction of minimal curvature u < of the free energy. The directions of the two orthogonal eigenvectors u < and u > define a set of orthogonal eigen-observables (ρ < , ρ > ) corresponding to a rotation of the set (ρ n , ρ p ) with an angle α = (u < , u n ):

ρ < = ρ n cosα + ρ p sinα ρ > = ρ n sinα -ρ p cosα (21)
At each point of the system evolution along the axis ρ < , the inverse curvature along the orthogonal axis measures the fluctuation of ρ > . Since the direction of phase separation u < is essentially isoscalar (α ≈ π/4), the observable ρ < is close to the total system density while

ρ > ≈ ρn-ρp √ 2 = ρ 3 √
2 approximately gives the isospin content of the system. In that sense, the fluctuations of ρ > are related to the isospin fluctuations of the system. The actual isospin fluctuations can be calculated projecting the fluctuations of ρ > on the ρ 3 direction.

If the process occurs at equilibrium, phase separation follows Gibbs construction rules, which constrain the two phases on the coexistence curve. Thus, the physical equilibrium fluctuations for the system are associated with the curvature of the free-energy surface f eq (ρ n , ρ p ) in the direction ρ coL (ρ coG ), given by the tangent to the coexistence curve at the liquid (gas) point.

Since phase separation is a simple straight line, these final fluctuations are nothing but the extrapolation to the coexistence border of the f eq inverse curvature in the u eq > direction, calculated at the initial point (ρ n , ρ p ) . The different curvatures associated with the two phases in coexistence are represented in fig. 6 as a function of the isovector chemical potential µ 3 = µ n -µ p , in order to have a common abscissa for liquid and gas. The projections in the ρ 3 direction are given in the lower part of the figure. The curvature of f eq is almost constant inside the coexistence zone (dashed lines). Because of the characteristic shape of the coexistence curve (cf fig. 3) this curvature is reduced (increased) at the liquid (gas) border. It is represented for both phases in full lines. This very general result is only due to the free-energy structure of the two phases, and we may expect it to stay qualitatively valid independently of finite-size effects and phase-separation mechanism. The extreme values of µ 3 correspond to the critical points, which are defined by a zero curvature along the coexistence curve. At these points, both curvatures in ρ coL and ρ coG become zero, fourth order terms have to be considered in the expansion (18), and the gaussian approximation (17) breaks down. Since the coexistence curve is closed, there is a value of asymmetry for which the coexistence line is tangent to the isoscalar direction, leading to a divergence of the ρ 3 -projected curvature on the liquid side.

However, this still happens close to the critical point, where the gaussian approximation is not enough.

Far from the critical point, these inverse curvatures can be associated with a fluctuation through eq.( 19). These fluctuations are largely independent of the system temperature inside the coexistence zone. They can be strictly interpreted as isospin fluctuations only for symmetric matter (µ 3 = 0), where the ρ coL and ρ coG directions are purely isovector. For symmetric matter, the liquid fluctuation appears much larger than the gas one. Because of the actual shape of the coexistence line, the isospin part of the liquid fluctuation is a rapidly decreasing function of the global asymmetry, as shown in the lower part of fig. 6. Conversely, the gas branch presents a fluctuation increasing with the asymmetry.

As far as physical fragmenting systems are concerned, fluctuations are affected by conservation laws which are not accounted for in this grand-canonical approach : therefore, the fluctuations represented in the lower part of fig. 6 can only be qualitatively interpreted as widths of fragment isotopic distributions.

Let us now focus on the difference between equilibrium fluctuations and those expected in spinodal decomposition. In both cases we will use the same physical picture so that results can be qualitatively compared.

In the preceeding discussion we have assumed that phase separation is entirely deter-mined by equilibrium rules. If, out of equilibrium, spinodal decomposition is the dominant mechanism, fluctuations in the two dimensional (ρ n , ρ p ) plane grow in time from the instability region until they cause the decomposition of the system. Then isospin fluctuations should be calculated as a dynamical variable given by the projection over the isovector axis of the two dimensional fluctuations varying in time. This approach is followed in numerical codes [5,26], as well as in simplified analytical models like in ref. [27], which follows the time evolution of unstable modes under the influence of a stochastic mean field in the linear approximation. Both approaches give similar predictions for the isotopic widths [27], and suggest narrower distributions than in experimental data and in the case of fragmentation in statistical equilibrium [26].

In these approaches, the initial condition is given by statistical equilibrium inside spinodal region of the observables which are not associated with the instability. The initial fluctuation of these observables is neglected. However, the inverse curvature along the axis orthogonal to the instability direction u 0 < is never zero, giving rise to an initial fluctuation of ρ > . We now want to estimate this effect. Consistently with the reasoning of section IV above, we assume that this fluctuation is mainly determined at the onset of phase separation, which is justified in ref. [28].

As we can see from fig. 7, this inverse curvature is not negligible but rather comparable to the case of equilibrium. In this figure, we compare the cases of equilibrium and spinodal decomposition for the three different Skyrme parameterizations. We consider systems along the axis joining the critical points. The abscissa is the isovector density of the global system, ρ 3 = ρ n -ρ p . For each point, we represent the ρ 3 projection of the different inverse curvatures.

For spinodal decomposition, the relevant free energy is f 0 , and its inverse curvature is calculated in u 0 > direction. At equilibrium, the concave intruder in f 0 is eliminated by Gibbs construction, and the inverse curvature of the resulting free energy f eq is computed in u eq > direction for the global system, which can be directly compared to the convexity of the unstable homogeneous system.

The main features obtained do not depend on Skyrme parameterization. We can see that, inside spinodal region, the convexification of the free energy does not modify drastically the curvature in the isospin direction. If spinodal decomposition is the dominant mechanism of phase separation, the isospin fluctuation tends to decrease with increasing global asymmetry, as in the equilibrium case. A complete dynamical treatment of the problem is needed to have quantitative predictions for the final fragment isotopic widths [26,27]. However in the minimal hypothesis of a completely diabatic dynamics that propagates the initial fluctuation over the free energy surface, we can expect the absolute isotopic width of fragments issuing from spinodal decomposition to be comparable to the case of coexisting phases at equilibrium. The high asymmetry region where this trend is inversed is not physical, since the gaussian approximation breaks down there for the liquid fraction in the equilibrium case.

To get more insight on this issue, we can compare the (initial) fluctuations expected for the spinodal mechanism with the predictions of the equilibrium SMM model in the grancanonical approximation. Within this formalism the production probability of an isotope composed of Z protons and N neutrons reads P β,µn,µp ∝ e β(-F β (N,Z)+µnN +µpZ)

Let us denote A 3 = V ρ 3 = N -Z. The isospin variance of the distribution ( 22) is

σ eq A 3 2 =< A 2 3 > -< A 3 > 2 = V 2 σ eq ρ 3 2
. Using the same saddle point approximation we used to get eq.( 19), this is related to the free-energy curvature by

σ eq A 3 2 = 1 β ∂ 2 F ∂A 2 3 -1 = A 2βC sym (A) (23)
where C sym (A) is the symmetry energy of the produced fragment. Denoting ρ eq the freezeout fragment density, we have C sym (A) = C sym (ρ eq ).

In the minimal hypothesis that the initial fluctuations σ sd A 3 2 i are not amplified (per unit volume) by the successive dynamics of phase separation :

σ sd A 3 2 f /V f = σ sd A 3 2 i /V i = V i σ sd ρ 3 2 i , (24)
the final variance in the case of spinodal decomposition σ sd A 3 2 f reads (cf eq.( 19) and ( 20)) :

σ sd A 3 2 f = V f ρ i 2βC sym (ρ i ) = A 2βC sym (ρ i ) ρ i ρ sd f (25)
where ρ i is the density at the onset of fragment formation and ρ sd f the freeze-out fragment density in a spinodal decomposition scenario. Equations ( 23) and ( 25) lead to the relation :

σ sd A 3 2 f σ eq A 3 2 =
C sym (ρ eq ) C sym (ρ i )

ρ i ρ sd f ( 26
)
Using the standard representation [29] C sym (ρ) = C sym (ρ 0 )(ρ/ρ 0 ) γ and considering ρ sd f ≈ ρ eq ≈ ρ 0 , we get

σ sd A 3 2 f σ eq A 3 2 ≈ ρ ρ 0 1-γ (27)
which can be even greater than 1 for an asy-stiff equation of state.

Considering that in the dynamics of spinodal decomposition isovector density fluctuations will also be (slightly) increased from the initial condition [27], eq.( 27) confirms that the fluctuations associated with the two scenarios are expected to be comparable.

VI. CONCLUSION

In this paper, we have presented isospin properties of the phases formed in nuclearmatter liquid-gas transition. The characteristics of phase separation are deduced from the free-energy curvature properties, f corresponding to the constrained entropy at a fixed temperature β. Under a critical temperature, the free energy of the homogeneous system f 0 presents a region of abnormal curvature, where phase separation is favorable. Phase equilibrium is obtained according to Gibbs construction, determining an equilibrated free energy f eq . Phase properties are then deduced from the free-energy curvature matrix, studying f eq for a system at equilibrium and f 0 for a system undergoing spinodal decomposition.

The direction of phase separation is given by the eigen-vector u < linked to the lower eigenvalue of the curvature matrix. In both equilibrium and spinodal decomposition, it leads to isospin fractionation, with a liquid more symmetric than the global system. In the direction orthogonal to phase separation, the curvature is linked to fluctuations affecting phase composition. This direction being close to ρ 3 direction, it is strongly linked to isospin fluctuations. Spinodal decomposition is predicted to give a stronger average fractionation and at least comparable or even higher isospin fluctuations with respect to a phase separation at equilibrium. Isoscaling observables as well as isotopic distributions are expected to be sensitive to such effects. If fragmentation can be associated with a spinodal decomposition, we qualitatively expect i) low apparent values of the symmetry energy coefficient extracted from isoscaling analyses, and ii) isotopic widths comparable or even larger than in the case of the statistical model. Finally we would like to stress that the present work concerns an idealized study of bulk matter, which is likely subject to significant modification when the nuclear surface is taken into account. The pursuit of this issue would be a very interesting topic for further study.

FIG. 1 :

 1 FIG. 1: Constrained entropy at T = 10M eV , corrected by βµ s ρ in order to emphasize curvature properties. It is shown on the left in perspective, and on the right as a contour plot. On the right, coexistence and spinodal boundaries are added as thick lines. Outer curve: coexistence. Straight lines link selected couples of phases at equilibrium. Inner curve: spinodal. Dots: critical points.

  The spinodal region is by construction inside coexistence. Coexistence and spinodal boundaries are tangent at the critical points. These curves are represented on the right part of fig.1. The phase diagram of the three different Skyrme forces is compared for the same temperature T = 10 MeV on the left part of fig.2. Sly230a and SGII give close results, while SIII shows an atypical behaviour, especially, as expected, on the asymmetric parts of the diagram.

FIG. 2 :

 2 FIG.2: Coexistence and spinodal curve in the density plane at T = 10M eV . Left: comparison between different Skyrme parameterizations : Sly230a (thick lines), SGII (thin lines) and SIII (dashed lines). Right: illustration of phase separation for a system of proton fraction Z/A = 0.3 inside instability region (cf section III).

FIG. 3 :

 3 FIG.3: Global representation of phase-separation directions through the density plane, for three different Skyrme parameterizations. Coexistence and spinodal boundaries are represented as in fig.2. A few selected equilibrium directions are given by straight lines joining the high and lowdensity branches of the coexistence curve (these straight lines are everywhere tangent to the equilibrium direction u eq < ). Phase-separation directions in spinodal decomposition are represented by a selection of lines which are everywhere tangent to the instability direction u 0 < . They are drawn as full lines inside the spinodal region where they correspond to an instability direction, and prolongated in dashed lines ouside the spinodal, where there is no more local instability.

FIG. 4 :

 4 FIG.4: Phase-separation directions δρ 3 /δρ at equilibrium (dash) and in spinodal decomposition (solid) as a function of the asymmetry ρ 3 /ρ for three different Skyrme parameterizations at different values of ρ : below, at and above the associated critical density. The instability direction is continued by the lowest curvature direction (dotted line). The dot-dash straight line gives the direction of constant Z/A while the double-dot-dash horizontal line gives the isoscalar direction δρ 3 = 0. Small dots give spinodal boundaries and large ones coexistence region boundaries.

FIG. 5 :

 5 FIG. 5: Influence of isospin fractionation on C sym measurement at T = 4M eV . Left side: phaseseparation directions. Top-left: separation direction at a typical density of fragment formation ρ = ρ ini = 0.05f m -3 , as a function of the global proton fraction for equilibrium (dashed line) and spinodal instability (full line). Vertical lines indicate the two systems S1∼ 112 Sn, S2∼ 124 Sn. Bottom-left: representation of the corresponding phase separations in the density plane. Right side: consequence on symmetry-energy measurement. Top-right : average isotopic composition of the liquid-like fraction as a function of fragment density at separation, for the two systems S1 and S2. Bottom-right: Symmetry energy coefficient C sym compared to the approximate expression obtained neglecting fractionation C 0 sym , as a function of separation density, for equilibrium (dashed line) and spinodal instability (full line). In the case of equilibrium, the liquid position is fixed in the density plane: it is indicated by the empty circles.

5 FIG. 6 :

 56 FIG. 6: Inverse curvature of the equilibrated free energy f eq as a function of the isovector chemical potential µ 3 . Dashed lines: curvature taken in the direction u eq > for points inside the coexistence region along the axis ρ = ρ c . Thick (thin) lines: curvature taken in the direction u coL (u coG ) tangent to the liquid (gas) branch of the coexistence curve. Left part: T = 4M eV . Right part: T = 10M eV . Lower part: projection on ρ 3 axis of the different curvatures considered.

FIG. 7 :

 7 FIG.7: Study of isospin fluctuations, compared at equilibrium and for spinodal decomposition for three different Skyrme parameterizations as a function of the total isovector density. The represented quantity is the inverse curvature of the free energy taken in the relevant direction and projected on ρ 3 axis. Global systems are taken along the axis ρ = ρ c . The system undergoing spinodal decomposition is represented by thin lines (curvature of f 0 taken in u 0 > direction). The curvature in u eq > direction is shown by the dashed lines. Dots give the position of critical points.


Nicolas P Rougier

: Glumpy

Fast

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

Glumpy is a python library for scientific visualization that is both fast, scalable and beautiful. Glumpy offers an intuitive interface between numpy and modern OpenGL. Glumpy is a laboratory for experimenting ideas to be integrated into the vispy project that lives at vispy.org, But glumpy can also be used as a standalone library.

  Numpy arrays are tracked automatically such that any change on the CPU is reflected automatically on the GPU, making the integration seamless. glumpy provides more than 100 examples showing how to achieve this or that visualization, such as antialiasing, text rendering, isolines computation, image filtering, cartographic projection, realtime signals, FLUID simulation, SVG rendering, etc. antialiased lines and text... spherical harmonics on the GPU text is rendered using signed distant fields did you know that a set of cones seen from above is a voronoi diagram ? import numpy as np from glumpy import app, collection window = app.Window(1024,1024, color=(1,1,1,1)) @window.event def on_draw(dt): window.clear() points.draw() points = Collection.points() points.append(np.random.normal(0.0,0.5,(1000000,3))) app.run() ...with antigrain quality nicolas.Rougier@inria.fr there's more at glumpy.github.io