Azaroual M Brgm
email: g.picot@brgm.fr

NUMERICAL MODELING OF FLUID-ROCK CHEMICAL INTERACTIONS DURING CO 2 SATURATED WATER INJECTION INTO A SANDSTONE RESERVOIR, USING THE MARTHE-REACT CODE

This study describes the numerical modeling of fluidrock chemical interactions during CO 2 -saturated water injection into a sandstone reservoir, using the MARTHE-REACT code (Thiéry et al., 2009, presented in the same TOUGH symposium). The reservoir is modeled in a 2D radial geometry with a thickness of 100 m, divided into three horizontal layers. The layers are characterized by contrasting hydrodynamic properties (porosity between 4 to 7%, permeability between 2×10 -16 m 2 and 2.3×10 -15 m 2 ), but all layers initially contain the same material (mineralogical assemblage and water composition). CO 2 -saturated water is injected vertically over the entire thickness of the reservoir.

The simulation results confirm the high reactivity of CO 2 -saturated water with respect to carbonate and sulphate minerals, which are totally dissolved in the water-saturated zone of the near-well field. Reactions occur to a larger extent in the upper layer, where the permeability is higher. Dissolution of calcite, dolomite, and anhydrite induces some pH buffering and a significant increase in porosity, reaching 41% in the first few meters of the water-saturated zone after ten years. Just behind the dissolution front, one can observe precipitation of dolomite, kaolinite, and quartz. Albite dissolution presents the first step in destabilization of the alumino-silicates. This can be explained by the slow kinetics of alumino-silicate mineral dissolution compared to the faster kinetics of carbonates. Consequently, the porosity evolution of the reservoir is essentially constrained by the magnitude of carbonate-and sulphate-mineral dissolution-precipitation reactions.

This study is a first step in the assessment of CO 2saturated water leakage from the reservoir to the surface in a sedimentary basin, where physicalfunctionality approaches implemented in the MARTHE simulator are essential for modeling complex hydrodynamic systems (supporting the link between hydrogeology and hydrology, including, for instance, mass exchanges with the atmosphere and river subsystems).

INTRODUCTION

We describe physical and chemical phenomena occurring during the injection of acidic solutions (CO 2 -rich solutions) with the establishment of reactive and mobile fronts.

Previous studies of geochemical processes occurring during CO 2 injection in an aquifer forecast the existence of five contrasting reactive zones extending radially from the injector (Figure 1): a full drying-out of the medium (Zone 5); a highly saline aqueous solution (Zone 4); a multiphase system with pH controlled by CO 2 (Zone 3); a saturated and acidified (CO 2 -rich solution) aqueous solution (Zone 2); and an aqueous solution at equilibrium with reservoir minerals at reservoir pH (Zone 1) (Gaus et al., 2008;André et al., 2007;Azaroual et al., 2007).

Figure 1. Conceptual diagram of the geochemical

processes believed to occur during CO 2 injection in an aquifer. r 1 , r 2 , r 3 , r 4 are the different radii; S G is gas saturation and S L is liquid saturation (from Gaus et al., 2008;André et al, 2007 andAzaroual et al, 2007).

These reactions may change the porosity and connectivity of the porous/fractured medium, and therefore have an impact on the performance of the injection. This study is focused on Zone 1 and Zone 2, which have single-phase flow. In our study, we assume the point zero is located at the interface between supercritical CO 2 and water (Zone 3-Zone 2). We analyze the main physical chemical processes involved with injection of CO 2 and saturated water into the deep Keuper aquifer in the Paris Basin (France), using a reactive transport modeling approach.

NUMERICAL TOOLS

The simulations presented in this study were performed using the MARTHE-REACT code (Thiéry et al., 2009, presented in the same TOUGH symposium), which was formed by coupling the MARTHE code (flow and transport in porous media) with the chemical part of the TOUGHREACT code.

TOUGHREACT reactive code

TOUGHREACT is a nonisothermal multicomponent reactive fluid flow and geochemical transport simulator (Xu et al, 2004) that can be applied under various thermohydrological and geochemical conditions of pressure (1 bar to several hundred bars), temperature (0-300°C), water saturation (dry to saturated), and ionic strength (up to 6 mol/kg H2O for NaCl-dominant solutions). Note that only the chemical part of TOUGHREACT is taken into account for the coupling with MARTHE.

MARTHE transport code

The MARTHE code (Modelling Aquifers with Rectangular cells, Transport and Hydrodynamics) is designed by BRGM for the hydrodynamic and hydrodispersive modeling of groundwater flow in porous media (Thiéry, 1990). Hydrodynamic calculations are carried out by a finite volume algorithm based on a matrix resolution of conjugate gradients with Choleski preconditioning. Several methods can be used for transport calculation (total variation diminishing, finite volumes, MOC, or random walk).

MARTHE-REACT coupling approach

MARTHE has already been coupled to specific chemical simulators (SCS) (Thiéry, 1995). The approach used to couple it with the TOUGHREACT chemical simulator is essentially the same as before.

Flow calculations and mass and energy transport are performed by the MARTHE code with no modifications. The chemical simulations are performed by the chemical part of the TOUGHREACT code. The coupling sequence is the following at each time step: (1) coupled flow calculation in aquifers and river networks integrating evaporation, infiltration, runoff, overflow, etc.;

(2) transport in aquifer and river networks of every primary dissolved chemical species; (3) chemical reactions in each aquifer cell and river reach using the concentrations at the end of the time step. There is no iteration; hence, it is a SNIA (Sequential Non Iterative Algorithm) scheme. The important point is that with this sequential scheme, after transport, the chemical-reaction calculations are independent in each cell.

INJECTION OF CO 2 -SATURATED WATER INTO THE KEUPER AQUIFER

Keuper aquifer properties

The study site is located in the Paris Sedimentary Basin, which includes deposits from the Mesozoic to the Cenozoic on a pre-Permian basement. The Keuper formation within the Paris Basin was recently studied for its CO 2 storage feasibility (Bonijoly et al., 2003;Delmas, 2007). Several aquifers in the east and south-central parts of the Paris Basin could be affected by CO 2 storage; we focus our study on the Donnemarie sandstone formation, a multilayered reservoir located at a depth of around 2,000 m with a thickness of 100 m. The reservoir, which is entirely sealed by anhydritic clays, is still exploited for hydrocarbons deposits and natural gas storage. Three layers (an upper, middle, and lower part) can be distinguished, with thicknesses of 35 m, 20 m, and 45 m, respectively. The upper part of the reservoir is the most permeable layer, whereas the porosity of the three parts is similarly low (Table 1). Water salinity is in the range 30-180 g kg -1 H2O . Temperature and pressure conditions in the reservoir (1,500 to 3,000 m depth) vary respectively from 70 to 120°C and 200 to 300 bars (Bonijoly et al., 2003).

Geometry

The Keuper aquifer (Donnemarie) is modeled as a sandstone reservoir represented by a cylindrical geometry with three layers, here representing the upper, middle, and lower parts of the reservoir (Figure 2). The grid is composed of 300 cells: along the Y axis, 1 cell with a radius of 0.2 m contains the vertical injection well, whereas 99 cells represent a radius following a logarithmic progression up to 100 km; the Z axes is divided by 3 cells corresponding to the upper (35 m), middle (20 m), and lower (45 m) parts of the Donnemarie sandstone. Each layer is assumed to be homogeneous, isotropic, and without slope.

Basin mineralogy

Formation water

The composition of the Keuper aquifer formation water is calculated at equilibrium with the mineralogy cited above. For the following simulation, the reservoir is assumed to be saturated by saline brine (≈ 150 g/kg H2O ), with the chemical composition given in

3 2 2 2 CO H O H CO       3 3 2 HCO H CO H
The resulting CO 2 -saturated solution, with a concentration about 1.3 mol CO 2 kg -1 H2O and a pH of 3.15, is injected into the reservoir (Table 2). A tracer is also injected in the reservoir to show the velocity of a nonreactive compound in the different layers (its concentration is 5×10 -5 mol.kg -1 H2O ).

   n n n n 1 A k r   
with k n the rate constant (mol/m 2 /s), A n the specific reactive surface area (m 2 ),  n the mineral saturation ratio (-) and  and  power terms defined from experiments. Calcite and iron hydroxides are assumed to dissolve-precipitate at equilibrium, because their reaction rates are quite rapid. Other minerals are kinetically controlled by acid, neutral, alkaline (1), and carbonation (2), mechanisms depending on temperature (Lasaga et al., 1994;Palandri and Kharaka, 2004):

H n H H a H 25 nu a nu 25 a 298.15 1 T 1 R E exp k 298.15 1 T 1 R E exp k k                                                                   ) 2 a 298.15 1 T 1 R E exp k ) 1 a 298.15 1 T 1 R E exp k OH aq 2, OH n CO OH a OH 25 n H OH a OH 25
with k 25 the rate constant at 25°C, E a the activation energy, R the gas constant, T the absolute temperature, a the activity of species, and n a power term, where the subscripts nu and H indicate neutral and acid mechanisms respectively. The acid catalysis is considered only for mineral dissolution. Table 3 provides all parameters defining the kinetic rate constrained for the simulation. For several minerals, the volume fraction is null, because minerals are not initially present in the reservoir, but allowed to precipitate.

Simulation scenario

The scenario is a hypothetical injection of a CO 2saturated water conducted to test the reactive geochemical model in Zone 2 (acidified zone) and observe the dissolution-precipitation phenomena on Keuper mineralogy. Injection rate of CO 2 -saturated water is chosen to have an increase in pressure lower than 20% of the initial pressure (200 bars) in each layer. Deduced injection rates correspond to 8.4×10 -4 m 3 /s, 6×10 -5 m 3 /s, and 9.45×10 -5 m 3 /s into upper, middle, and lower layers respectively. An amount of 300,000 m 3 of CO 2 -saturated water flows over 10 years, which corresponds to 18,720 tons of dissolved CO 2 considering 1.3 mol CO 2 kg -1 H2O . A hydrostatic pressure is considered initially and is fixed in the outermost column of the grid during the simulation. No exchanges of fluids are considered under and above the system.

Simulation results

Hydrodynamics

The tracer included within the CO 2 -saturated water is represented on a profile (Figure 3) from 0 to 350 m after 10 years of injection. Its behavior confirms different flow velocities for the three layers, with the fastest transport in the upper layer (the most permeable). After 10 years, originally injected CO 2saturated water reaches a radius of 150 m in the upper layer, 70 m in the middle layer, and 45 m in the lower layer. We can expect preferential geochemical reactivity in these zones.

Geochemical reactivity

As original injected CO 2 -saturated water reaches a radius of 150 m in the upper layer, we could expect a pH of 3.15 in this area. But the pH of 3.15 exists only in the first few meters around the injection well (Figure 4). Further into the upper layer (radius 3-150 m), the pH is 4.8, which indicates that geochemical reactivity is buffering the pH. Estimated changes in porosity and permeability are located in the same area as well (Figure 5). Figure 5. Profile of potential impact on porosity and permeability (red zone) after 10 years of CO 2 -saturated water injection. Maximal porosity is estimated around 0.41 and maximal permeability is deduced equal to 7×10 -13 m 2 .

For a better understanding, we focus on the geochemical reactivity of the upper layer (Figure 6). As observed previously, the pH of 3.15 extends a few meters (0 to 2.5 m), forming an acidic front: Zone 2a (A). The total dissolution of dolomite, calcite (B), and anhydrite (C) is observed in the same area.

                    2 2 3 2 4 2 3 Mg Ca HCO 2 H 2 Dolomite Ca HSO H Anhydrite Ca HCO H Calcite
Dissolution of albite (D) is observed in lower quantities as a result of smaller kinetic rates.

        Na Al SiO H 3 O H 4 H 4 Albite 3 4 4 2
The potential impact of mineral dissolution on porosity is estimated (A), but not taken into account in the flow calculation. A maximal increase in porosity (41% from 6% initially) is estimated, which corresponds to the total dissolution of carbonate and sulphate minerals.

 

Precipitations of kaolinite and quartz have a minor impact on porosity changes compared to dissolution of carbonates. But at this time scale (10 years), porosity changes are limited by the amount of carbonate, because of the slow reactivity of silicates. Further, in the buffered zone (from 150 m to 100 km), pH increases to reach the pH of formation water, and any dissolution-precipitations are observed because of the equilibrium between minerals and formation water (zone 1 = initial zone). The simulation results confirm the high reactivity of CO 2 -saturated water with respect to the carbonate minerals (André et al., 2007), which are totally dissolved in the water-saturated zone of the near-well field. Reactions occur on a larger extent in the superior layer, where the permeability is higher. Dissolution of calcite, dolomite, and anhydrite induces some pH buffering and a significant increase in porosity, reaching 41% in the first meters of the water-saturated zone after ten years, when the dissolution is total. Albite dissolution presents the first step in destabilization of the alumino-silicates in lower quantities. This can be explained by the slow kinetics of dissolution for alumino-silicate minerals compared to that for carbonates. This is in accordance with observations by Audigane et al. (2007), who simulate CO 2 injection into sandstone. Just behind the acidic zone, one can observe precipitation of dolomite, kaolinite, and quartz. The porosity evolution of the reservoir is essentially constrained by the magnitude of carbonate and sulphate mineral dissolution, and the impact on injectivity still depends on the location of mineral dissolution-precipitation (connected pores or not).

CONCLUSION

A 2D-radial numerical model of CO 2 -saturated water injection into the sandstone Keuper aquifer in the Paris Basin has been presented using the MARTHE-REACT code. The simulation was performed over 10 years of continuous injection considered as a high reactivity level compared to supercritical CO 2 injection as studied by André et al (2007). Three reactive zones were identified: 1) An acidic zone where dolomite, calcite, and anhydrite are totally dissolved, and where albite starts to dissolve from the acidic conditions (pH Acidic Zone

Buffering Zone

Initial Zone 1 2b 2a = 3.15), with a maximal porosity increase of 41%, 2) A buffered zone where carbonate and sulphate minerals buffer the acidic injection (pH = 4.8) and where transported ions allow precipitation of kaolinite, quartz, and dolomite with minor impact on porosity, 3) An initial zone where minerals are in equilibrium with formation water (pH = 6.46) without changes in porosity.

The results obtained confirm similar findings from other studies on the reactivity of carbonate and sandstone aquifers, and provides further confidence in the coupling tool.

This study is a first step before the assessment of CO 2 -saturated water leakage from a reservoir to the surface in a sedimentary basin, where functionalities of the MARTHE simulator are essential for modeling complex hydrodynamic systems (providing a link between hydrogeology and hydrology, including, for example, atmosphere and river subsystems).

Figure 2 .

 2 Figure 2. Two-dimensional radial geometry of the Keuper reservoir model with three layers representing upper, middle and lower parts of Donnemarie sandstones 0 100 km

Figure 3 .Figure 4 .

 34 Figure 3. Profile of tracer after 10 years of CO 2 -saturated water injection. Initial concentration of tracer is equal to 5 10 -8 mol.kg -1 H2O and injected tracer concentration is equal to 5 10 -5 mol.kg -1 H2O .

Figure 6

 6 Figure 6 Variation of porosity, pH, and mineral volume fractions (Vmin/Vs) in the upper layer of Donnemarie aquifer from the injection cell to the outermost after 10 years of CO 2 -saturated water injection. Curves are plotted on semilog scales.

Table 1 .

 1 Physical properties of the Keuper aquifer (Donnemarie).

			Reservoir part
		Upper	Middle	Lower
	Porosity (-)	0.06	0.04	0.07
	Permeability (m²)	2.3 10 -15	0.3 10 -15	0.2 10 -15

Table 2

 2

	at a temperature of
	100°C and 200 bars pressure.
	CO 2 -saturated injected water
	CO 2 -saturated injected water is created from the
	equilibrium between initial dilute formation water
	with a pure CO 2 phase at 100°C and a pressure of 200
	bars (equivalent to a CO 2 fugacity of 123.1 bars). As
	mentioned in André et al. (2007), this reaction
	produces bicarbonate ions, which gives an acid
	solution (pH of 3.15).

Table 2 .

 2 Chemical composition of water from the Keuper aquifer (formation water) and of CO 2 -saturated water injected into the Keuper aquifer (mol.kg -1 H2O )

		Formation	CO2-saturated
		Water	Water
	Temperature	100°C	100°C
	pH	6.46	3.15
	C	3.343e-3	1.3069
	Ca	3.068e-2	3.071e-5
	Mg	8.064e-3	8.061e-6
	Na	2.468	2.468e-3
	K	1.376e-3	1.376e-6
	Fe	1.844e-13	1.010e-13
	H4SiO4	8.93e-4	8.932e-7
	SO4	6.183e-2	8.181e-5
	Al	3.028e-7	3.028e-10
	Cl	2.42	2.42e-3
	Tracer	5.063e-8	5.00e-5
	Geochemical approach
	Kinetic rates are taken into account to examine the
	different mineral behaviors. Mineral (n) dissolution
	(+) -precipitation (-) rate (r

n ) in mol/s is expressed by

Lasaga et al. (1994)

:

Table 3 .

 3 Kinetic rate parameters of each mineral for neutral and acid mechanism and associated volume fraction *Dissolve and precipitate at equilibrium.

				Acid		Neutral		Carbonation
	Minerals volume fraction	Mechanism		Mechanism		Mechanism
		(Vmin/Vs)	Log k	Ea	n	Log k	Ea	Log K	Ea	m
	Dolomite Siderite	0.25 0.0	-3.19 36.1	0.5	-7.53	52.2	-5.11	34.8	0.5
	Dawsonite	0.0	-	-	-	-7	62.8	-	-	-
	Magnesite	0.0	-6.38 14.4	1	-9.34	23.5	-5.22	62.8	1
	Calcite*	0.09							Alkaline
					-				mechanism
	Fe(OH)3*	0.009
								Log k	Ea	n
	Albite	0.08	-10.16	65 0.457 -12.56	69.8	-15.6	71	-0.572
	Kaolinite	0.0025	-11.31 65.9 0.777 -13.18	22.2	-17.05	17.9	-0.472
	Illite-Mg	0.01	-10.98 23.6 0.34 -12.78	35.0	-16.52	58.9	-0.4
	Quartz	0.5285				-13.4	90.9
	Halite	0.0	-	-	-	-0.21	7.4	-	-	-
	Anhydrite	0.03	-	-	-	-3.19	14.3

  SiO 4 in solution, allowing precipitation of kaolinite and dolomite just after the acidic zone (2.5 to 3 m).

	These dissolutions release Ca 2+ , Mg 2+ , HCO 3 -, HSO 4 -, Al 3+ , and H 4 Kaolinite SiO H 2 O H Al 2 4 4 2 3    
	Further (3 to 150 m), we can observe a second area
	where the pH is buffered by carbonate and sulphate,
	and where kaolinite and quartz precipitate (zone 2b =
	buffering zone).
	H	SiO	Quartz	H	O
	4	4		2

ACKNOWLEDGMENT

This study is carried out in the framework of the PROCHE PUITS project co-funded by the French National Research Agency (ANR, in French).

André, L., P. Audigane, M. Azaroual, A. Menjoz, Numerical modeling of fluid-rock chemical interactions at the supercritical CO2-liquid interface during CO2 injection into a carbonate reservoir, the Dogger aquifer (Paris Basin, France), Energy Conversion & Management, 48, 1782-1797, 2007. Audigane, P., I. Gaus, I. Czernichowski-Lauriol, K. Pruess and T. Xu, Two-dimensional reactive transport modelling of CO2 injection in a saline aquifer at the Sleipner site, North Sea, American Journal of Science, 307, 974-1008, 2007. Azaroual, M., K. Pruess, C. Fouillac, Feasibility of using supercritical CO 2 as heat transmission fluid in the EGS integrating the carbon storage constraints, Engine Workshop 2, 1-4 April 2007, Volterra, Italy, p. 47, 2007. Azaroual, M., C. Fouillac, J.M. Matray, Solubility of silica polymorphs in electrolyte solutions, II. Activity of aqueous silica and solid silica polymorphs in deep solutions from the sedimentary Paris Basin, Chemical Geology, 140, 167-179, 1997. Bonijoly, D., J. Barbier, J.M. Matray, C. Robelin, C. Kervevan, D. Thiéry, A. Menjoz, C. Cotiche, B. Herbrich, Feasibility of CO 2 storage in geothermal reservoirs. Example of the Paris Basin, France, Report BRGM/RP-52349-FR, 2003. Delmas, J., Etude pétrophysique des aquifères gréseux triasiques dans le secteur d'étude de PICOREF (Sud Champagne, Bassin de Paris), Porosité et perméabilité mesurées sur carottes, Rapport d'avancement Projet PICOREF: Volet 3, 2007, restricted. El-ghali, M.A.K., S. Morad, H. Mansurbeg, M.A. Caja, M. Sirat, N. Ogle, Diagenetic alterations related to marine transgression and regression in fluvial and shallow marine sandstones of the Triassic Buntsandstein and Keuper sequence, the Paris Basin, France, Marine and Petroleum Geology, 26, 289-309, 2009. Gaus, I., P. Audigane, L. André, J. Lions, N. Jacquemet, P. Durst, I. Czernichowski-Lauriol, M. Azaroual, Geochemical and solute transport modelling for CO 2 storage, what to expect from it?, International journal of greenhouse gas control, 2, 605 -625, 2008. Lasaga, A. C., J. M. Soler, , J. Ganor,, T. E. Burch, K. L. Nagy, Chemical weathering rate laws and global geochemical cycles, Geochim. Cosmochim. Acta, 58, 2361-2386, 1994. Palandri, J.L., and Y.K. Kharaka, A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modelling, U.S. Geological Survey Open file report (of 2004-1068), 2004. Parkhurst, D.L., and C.A.J. Appelo, User's Guide to PHREEQC (Version 2)--A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations, Water Resources investigations report, 99-4259, 1999. Renoux, P. Aquifères : Etude diagraphique du puits Limoreau 101. Projet GEOCARBONE-PICOREF-ANR CO2 Edition 2005, 150-151, 2005. Thiéry, D., N. Jacquemet, G. Picot-Colbeaux, C. Kervévan, L. André and M. Azaroual, Validation of MARTHE-REACT coupled surface and groundwater reactive transport code for modeling hydro systems, TOUGH Symposium 2009, 14-16 September, submitted. Thiéry, D, Modélisation des écoulements avec interactions chimiques avec le logiciel MARTHE Release 5.5, Report BRGM/R38463-HYT/DR-95, 1995. Thiéry, D. Software MARTHE, Modelling of aquifers with a rectangular grid in transient state for hydrodynamic calculations of heads and flows, Release 4.3, Report BRGM/R32548, 1990. Xu, T., E. Sonnenthal, N. Spycher and K. Pruess, TOUGHREACT User's Guide: A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated

geologic media, Report LBNL-55460, Lawrence Berkeley National Laboratory, 2004.


Christophe Chatelain
email: christophe.chatelain@univ-lorraine.fr

DMRG study of the Berezinskii-Kosterlitz-Thouless transitions of the 2D five-state clock model

Keywords: numbers: 05.20.-y, 05.50.+q, 75.10.Hk, 05.70.Jk, 05.10.-a σ

The two Berezinskii-Kosterlitz-Thouless phase transitions of the twodimensional 5-state clock model are studied on infinite strips using the DMRG algorithm. Because of the open boundary conditions, the helicity modulus Υ 2 is computed by imposing twisted magnetic fields at the two boundaries. Its scaling behavior is in good agreement with the existence of essential singularities with σ = 1/2 at the two transitions. The predicted universal values of Υ 2 are shown to be reached in the thermodynamic limit. The fourth-order helicity modulus is observed to display a dip at the high-temperature BKT transition, like the XY model, and shown to take a new universal value at the low-temperature one. Finally, the scaling behavior of magnetization at the low-temperature transition is compatible with η = 1/4.

Introduction

The critical behavior of the two-dimensional XY model has been the subject of a huge literature. Because the symmetry group of the XY model is continuous, the Mermim-Wagner-Hohenberg theorem states that long-range order is destroyed by massless spinwave excitations in dimension d = 2 [START_REF] Mermin | Absence of ferromagnetism or antiferromagnetism in one-or two-dimensional isotropic heisenberg models[END_REF][START_REF] Hohenberg | Existence of long-range order in one and two dimensions[END_REF]. Therefore, in contrast to the Ising or 3 or 4-state Potts models, there cannot exist any ferromagnetic phase, where the symmetry would be spontaneously broken. However, at low temperature the two-dimensional XY model presents a critical phase that terminates at a topological Berezinskii-Kosterlitz-Thouless (BKT) phase transition where the coherence of spin waves is lost because of the proliferation of free topological defects [3,[START_REF] Kosterlitz | Ordering, metastability and phase transitions in twodimensional systems[END_REF][START_REF] Minnhagen | The two-dimensional coulomb gas, vortex unbinding, and superfluidsuperconducting films[END_REF]. BKT phase transitions were found in a variety of context, the superfluid transition of helium being certainly the most famous example of such a transition.

However, there exist situations where the Mermim-Wagner-Hohenberg theorem can be bypassed. For instance, the theorem extends to any kind of pair interaction but not to hard-core potentials [START_REF] Mermin | Crystalline order in two dimensions[END_REF]. As a consequence, a two-dimensional gaz of hard spheres crystallizes at low enough temperature. As the temperature is increased, the system undergoes two successive phase transitions separating, first the solid phase and an intermediate hexatic phase, and then, the hexatic and liquid phases [START_REF] Halperin | Theory of two-dimensional melting[END_REF]. The solid-hexatic phase transition is a BKT transition while the nature of the hexaticliquid transition is controversial: it was long believed to be a BKT transition too but recent Monte Carlo simulations provided evidence of a first-order phase transition [START_REF] Bernard | Two-step melting in two dimensions: First-order liquid-hexatic transition[END_REF].

Another way to escape from the Mermim-Wagner-Hohenberg theorem is to replace the U (1) symmetry of the XY model by a finite symmetry. This is the case of the q-state clock model, originally introduced by Potts [START_REF] Potts | Some generalized order-disorder transformations[END_REF]. The Hamiltonian is defined as

H 0 = - (i,j) cos 2π q (σ i -σ j ) , σ i ∈ {0, . . . , q -1} (1)
where the sum extends over all pairs of nearest neighbors of the lattice. It is invariant under the Z q group of cyclic permutations of the q states. The XY model is recovered in the limit q → +∞. For q ≤ 4, the two-dimensional clock model undergoes a unique second-order phase transition, like the Potts model. For q ≥ 5, the phase diagram shows, between the ferromagnetic and paramagnetic phases, an intermediate critical phase with quasi long-range order [START_REF] José | Renormalization, vortices, and symmetry-breaking perturbations in the two-dimensional planar model[END_REF]. These three phases are separated by two BKT phase transitions. The transition between the critical and paramagnetic phases occurs at a temperature roughly independent of q. The transition temperature of the second transition decays to zero as the number of states q goes to infinity, in agreement with the absence of ferromagnetic order predicted by the Mermim-Wagner-Hohenberg theorem for the 2D XY model. The situation is however not completely clear in the case q = 5. Early Monte Carlo simulations provided evidence of the existence of three distinct phases [START_REF] Tobochnik | Properties of the q-state clock model for q = 4, 5, and 6[END_REF]. The average complex magnetization was shown to display the expected Z 5 -symmetry in the ferromagnetic phase and a continuous U (1) symmetry in the critical phase [START_REF] Borisenko | Numerical study of the phase transitions in the two-dimensional z(5) vector model[END_REF]. If the low-temperature transition was recognized as a BKT transition, the question of the nature of the high-temperature transition was subject to controversy in the last years. The claim that the high-temperature transition could not be of BKT-type for q ≤ 6 [START_REF] Lapilli | Universality away from critical points in two-dimensional phase transitions[END_REF] was supported by a study of Fisher zeroes [START_REF] Hwang | Six-state clock model on the square lattice: Fisher zero approach with wanglandau sampling[END_REF]. However, recent Monte Carlo studies showed that the helicity modulus Υ 2 , defined as the response of the system to a twist [START_REF] Fisher | Helicity modulus, superfluidity, and scaling in isotropic systems[END_REF], displays a jump, as expected in the case of a BKT transition [START_REF] Baek | Non-kosterlitz-thouless transitions for the q-state clock models[END_REF][START_REF] Borisenko | Phase transitions in two-dimensional z(n) vector models for n > 4[END_REF][START_REF] Baek | Residual discrete symmetry of the five-state clock model[END_REF][START_REF] Kumano | Response to a twist in systems with Zp symmetry: The two-dimensional p-state clock model[END_REF]. The observation of a non-vanishing helicity modulus in the high-temperature phase was explained as the result of an inappropriate definition of Υ 2 in the case of the clock model [START_REF] Kumano | Response to a twist in systems with Zp symmetry: The two-dimensional p-state clock model[END_REF].

The helicity modulus is central in showing the BKT nature of the transition. However, because it is defined as a derivative of the free energy with respect to a twist, its computation is not straightforward in a Monte Carlo simulation. In the references mentionned above, the boundary-flip Monte Carlo method was employed [START_REF] Gottlob | The xy model and the three-state antiferromagnetic potts model in three dimensions: Critical properties from fluctuating boundary conditions[END_REF]. In this paper, we present new results for the second and fourth-order helicity modulus of the 5 and 6-state clock models obtained using the Density Matrix Renormalization Group (DMRG) algorithm [START_REF] White | Density matrix formulation for quantum renormalization groups[END_REF][START_REF] White | Density-matrix algorithms for quantum renormalization groups[END_REF][START_REF] Schollwöck | The density-matrix renormalization group[END_REF]]. Details about the numerical calculations and the definition of the helicity modulus is presented in the second section. The numerical data are analyzed in the third section. The scaling behaviour of the helicity modulus, of the fourth-order helicity modulus and of magnetization are compared to Kosterlitz-Thouless predictions. Conclusions follow.

DMRG algorithm, observables and convergence

In this section, the DMRG algorithm is very briefly presented. For more details about the implementation of this algorithm, the interested reader is invited to refer to the review [START_REF] Schollwöck | The density-matrix renormalization group[END_REF]. Our implementation and the simulation parameters are discussed. The definition of the estimator of the helicity modulus is given and finally, the convergence is studied.

Z M +1 (σ ′ 1 , . . . , σ ′ L ) = σ 1 ,...,σ L T (σ ′ 1 , . . . , σ ′ L ; σ 1 , . . . , σ L )Z M (σ 1 , . . . , σ L )
suggests to use a matricial notation: |Z M +1 = T |Z M where T is the transfer matrix acting in the q L -dimensional vector space of the spin configurations on the upper row of the two-dimensional system. On the right, schematic representation of the action of the effective transfer matrix T eff constructed by the DMRG algorithm. Each vector of the vector space on which acts T eff corresponds to a configuration of the two central spins and to effective degrees of freedom on both sides (represented as shaded rectangles).

From the Hamiltonian (1) on a square lattice, a transfer matrix T is easily defined in the vector space of spin configurations on the last row of a strip of width L (see figure 1). Because the hamiltonian couples only nearest-neighbouring spins, the transfer matrix can be expressed as a product of matrices, whose elements are the Boltzmann weights T (σ ′ , σ) = e βJ cos 2π q (σ ′ -σ) of each new single bond betwen the σ and σ ′ . This property is important for DMRG. The free energy of the clock-model on an infinite strip of width L is proportional to the logarithm of the largest eigenvalue of the transfer matrix T . However, since the spins take q possible values, the dimension of this vector space grows as q L , i.e. exponentially fast. In the case q = 5 we are interested in, this limits the strip widths that can be studied to small values, L = 9 or 10. In the Density Matrix Renormalization Group (DMRG) approach, also called Transfer Matrix Renormalization Group (TMRG) in the context of classical systems, the transfer matrix T is replaced by an effective transfer T eff matrix acting in a much smaller vector space. The vectors of the latter corresponds to the configurations of two spins surrounded by effective degrees of freedom, the left and right blocks. T eff can be constructed iteratively, i.e. by considering strips of increasing width L. For a small system, the eigenvector of the transfer matrix T associated to the largest eigenvalue can be computed exactly. Effective transfer matrices with a smaller dimensions are then constructed for the left and right blocks. As shown by White [START_REF] White | Density matrix formulation for quantum renormalization groups[END_REF][START_REF] White | Density-matrix algorithms for quantum renormalization groups[END_REF][START_REF] Schollwöck | The density-matrix renormalization group[END_REF], the best approximation, in the sense of mean-square deviation of the eigenvector, is to truncate the basis of the left and right blocks by keeping only the eigenvectors associated to the largest eigenvalues of the reduced density matrix of the block. The central spins are then absorbed into the left and right blocks and two new spins are usually introduced between them. The procedure is repeated until the desired strip width is reached. This algorithm, called infinite-size DMRG, is efficient for homogeneous system. In our case, we will introduce magnetic fields at the boundaries. Therefore, we proceeded differently: only one block, say left, is considered. The second block, say right, is reduced to a single spin. After the construction of an effective transfer matrix for the left block, new spins are not introduced between the former central spins but at the right of the right block. The procedure is iterated up to the desired strip width. The same algorithm is used to construct a sequence of effective transfer matrices for the right block. Once the final strip width L has been reached, the accuracy of the free energy can be further improved by applying several iterations of the finite-size DMRG algorithm. A sweep over the system is performed by increasing the length of one block and decreasing the length of the other one. We systematically performed six sweeps, in both directions. Even for the longest strips, the free energy was observed to reach a stable value (10 digits were considered), after three sweeps.

Since the DMRG algorithm gives directly access to the free energy, the helicity modulus Υ 2 can be calculated as

Υ 2 = L 2 ∂ 2 f ∂∆ 2 ∆=0 ( 2
)
where f (∆) is the free energy density of the system when a twist ∆ is introduced at its boundary. However, as emphasized by Kumano et al. [START_REF] Kumano | Response to a twist in systems with Zp symmetry: The two-dimensional p-state clock model[END_REF], spins are not continuous in the clock model but takes a finite number of states so the twist can only be a multiple of 2π/q. The usual definition of the helicity modulus should be replaced in this case by

Υ 2 = 2L 2 (2π/q) 2 f (2π/q) -f (0) . (3)
The twist is usually introduced in the system by replacing the interaction term in the Hamiltonian on one row by

-cos 2π q (σ i -σ j + 1) . (4)
With periodic boundary conditions, a spin wave is induced to accommodate this interaction term. Unfortunately, periodic boundary conditions are known to greatly deteriorate the convergence of the DMRG algorithm. We therefore introduced a twist by imposing magnetic fields at the two boundaries of an open system. The Hamiltonian now reads

H ∆ = - (i,j) cos 2π q (σ i -σ j ) - i∈SL cos 2π q σ i - i∈SR cos 2π q σ i -∆ (5)
where S L (resp. S R ) denotes the set of spins at the left (resp. right) boundary of the strip. Equivalently, the system can be seen as a strip of width L + 2 with spins frozen in the states σ i = 0 and σ i = q 2π ∆ on the first and last columns respectively. The helicity modulus is determined as the finite difference of the free energy (3). We also tested the same definition but with a larger twist ∆ = 4π/q. Following [START_REF] Minnhagen | Direct evidence of the discontinuous character of the kosterlitzthouless jump[END_REF], a fourth-order helicity modulus Υ 4 can be defined in the case of an infinitesimal twist ∆ from the Taylor expansion

f (∆) = f (0) + ∆ 2 2 Υ 2 + ∆ 4 4! Υ 4 + O(∆ 6 ) (6)
of the free energy density. We introduced the estimator

Υ 4 = 2 (2π/q) 4 f (4π/q) + 3f (0) -4f (2π/q) . (7)
Finally, since the Z q symmetry is broken, the spontaneous magnetization can be computed when ∆ = 0, i.e. when the magnetic fields at the two boundaries are aligned in the same direction.

In Monte Carlo simulations, error bars on any observable O are related to the fluctuations of O. According to the central limit theorem, the error on O is given by σ 2 O τ /N where σ 2 O is the variance of O, N the number of Monte Carlo steps and τ the autocorrelation time. Error bars are therefore purely statistical and a better accuracy is reached by increasing the number of Monte Carlo steps. Systematic deviations arise only when the system was not properly thermalized or when the dynamics is not ergodic. In DMRG calculations, there is no statistical fluctuation of free energy around its exact value but only systematic deviations due to the truncation of the vector spaces of the left and right blocks. These deviations are difficult to estimate. In the following, the convergence of the helicity modulus is studied at the inverse temperature β = 1/k B T = 1.08, i.e. T ≃ 0.925. This temperature roughly lies at the center of the critical phase of the 5-state clock model. It is therefore the most difficult case for numerical simulations. Since the correlation length is expected to diverge with the strip width L in this phase, the left and right blocks are strongly correlated through the two central spins. As a consequence, the reduced density matrix is observed to have a slowly decaying spectrum. Because the deviation of free energy is expected to be proportional to the sum of the eigenvalues of the reduced density matrix associated to the eigenvectors that have been discarded [START_REF] Legeza | Accuracy of the density-matrix renormalization-group method[END_REF], each truncation of the basis of the two blocks leaves an effective transfer matrix that may reproduce poorly these correlations if a too small number of states is kept.

On figure 2, the helicity modulus is plotted versus the number of states m kept after truncation. More precisely, m refers to the number of states describing both the left (resp. right) block and the left (resp. right) central spin after the renormalization step. The dimension of the full vector space is therefore m 2 . For the smallest strip widths L, the numerical data do not show any dependence on m. However, the free energy displays a small evolution at the 8th digit for L ≥ 64. Since the helicity modulus Eq. 3 is defined as the difference of two free energies, the deviation is larger for Υ 2 . Indeed, for L = 64, the free energies are of order O(10 4 ) while their difference is of order O(10 0 ). As a consequence, an evolution of Υ 2 with m is observed at the 4th digit for L = 64 and at the third one for L = 256. To quantify the systematic deviation of our data, we extrapolated them to the limit m → +∞. As observed for instance in the case of SU (N ) chains [START_REF] Fhringer | Dmrg studies of critical su(n) spin chains[END_REF], the free energy is well described by a power law decay f (m) = f (∞) + am -α with α ≃ 4. The helicity modulus Υ 2 is also reasonnably well fitted by such a law, though other laws, in particular Υ 2 (m) = Υ 2 (∞) + ae -m/m0 with m 0 ≃ 70, give sensibly similar results. At L = 256, the fit is morever hampered by the existence of oscillations that becomes more and more important at large strip widths. At L = 128 and β = 1.08, the helicity modulus can be extrapolated to the value 0.86450 (with both power-law and expoential fits) in the limit m → +∞. In the following, m = 325 states will be kept at each truncation. Since the helicity modulus takes the value 0.86453 for m = 325, the systematic deviation at L = 128 can be estimated to be of order O(10 -5 ). By the same procedure, the deviation is shown to be one order of magnitude smaller for L = 64. For L = 256, this systematic deviation is of order O(10 -4 ). In the following, numerical data for the 6-state clock will also be presented for comparison. For this model, m = 648 states were kept during the truncations but strip widths only up to L = 64 could be considered.

Numerical evidence of two BKT transitions

Helicity modulus

A clear signature of the BKT transitions is provided by the helicity modulus Υ 2 . The temperature dependence of Υ 2 is plotted on figure 3 for both the 5 and 6-state clock models. In the ferromagnetic phase, the twist imposed at the boundaries induces domain walls rather than spin waves. As a consequence, L f (∆) -f (0) is finite. Since the helicity modulus is proportional to L 2 f (∆) -f (0) , it diverges in the ferromagnetic phase, as can be observed on the figure. In the paramagnetic phase, the correlation length is finite so no (quasi) long-range order can propagate into the system and the difference L f (∆) -f (0) , as well as the helicity modulus, goes to zero as the strip width is increased. For large strip widths, the helicity modulus is perfectly equal to zero, in contradistinction to previous claims. In the critical phase, the helicity modulus is expected to take a size-independent value. If this behavior is clearly observed on figure 3 for the 6-state clock model, this is not yet the case for the 5-state model, even though the curves come closer to each other as the strip width is increased. The two expected jumps of helicity modulus are not observed yet but the tangent of the curve becomes more and more steep at temperatures close to the expected BKT transition temperatures T low BKT ≃ 0.90 and T high BKT ≃ 0.95. When the BKT transition point is approached from the non-critical phase, the helicity modulus takes a universal value which is predicted to be proportional to T BKT . In the case of the 2D clock model, these two universal quantities are [START_REF] Kumano | Response to a twist in systems with Zp symmetry: The two-dimensional p-state clock model[END_REF]:

Υ 2 (T high BKT ) = 2 π T high BKT , Υ 2 (T low BKT ) = q 2 8π T low BKT . (8)
These linear behaviors are plotted on figure 3. The different curves are expected to intersect each one of these two straight lines at the same point in the thermodynamic limit. As can be seen on the figure, and as already noted by Kumano et al. [START_REF] Kumano | Response to a twist in systems with Zp symmetry: The two-dimensional p-state clock model[END_REF], this is not yet the case in the range of strip widths that were considered. However, the accuracy of the DMRG data allows to determine the intersection of the helicity modulus with the two predictions [START_REF] Fhringer | Dmrg studies of critical su(n) spin chains[END_REF]. On figure 4, the temperatures at which occur these intersections are plotted versus 1/(ln L) 2 . Since the correlation length displays an essential singularity

0 5 10 15 -4 -3 -2 -1 0 1 2 (β -β high BKT ) (ln L/L 0 ) 2 (πβΥ 2 - 2) ln L/L 0 L = 16 L = 24 L = 32 L = 64 L = 128 L = 256
(πβΥ 2 - 2) ln L/L 0 (β low BKT -β) (ln L/L 0 ) 2 L = 16 L = 24 L = 32 L = 64 L = 128 L = 256
ξ ∼ e a|T -TBKT| -1/2 , T > T BKT (9)
the temperature shift T -T KT should indeed scale as 1/(ln L) 2 in the Finite-Size regime ξ ∼ L. The data are in good agreement with this statement. Linear fits lead then to the extrapolated BKT transition temperatures T low BKT ≃ 0.914 [START_REF] Hwang | Six-state clock model on the square lattice: Fisher zero approach with wanglandau sampling[END_REF] and T high BKT ≃ 0.945 [START_REF] Legeza | Accuracy of the density-matrix renormalization-group method[END_REF], in agreement with earlier estimates T low BKT ≃ 0.90514(9) and T high BKT ≃ 0.95147 [START_REF] Gottlob | The xy model and the three-state antiferromagnetic potts model in three dimensions: Critical properties from fluctuating boundary conditions[END_REF] [6] obtained by other, more accurate, techniques but not involving the helicity modulus.

A stronger prediction of Kosterlitz-Thouless theory is that the universal value of the helicity modulus is approached as

Υ 2 (T ) -Υ 2 (T BKT ) ∼ ln L L 0 -1 F Υ |T -T BKT |(ln L/L 0 ) 2 (10)
where F Υ is a universal scaling function. This prediction is tested on figure 5 for the two BKT transitions of the clock model. The collapse of the curves is reasonably good for a wide range of the scaling parameter. For comparison, the same curves are plotted for the 6-state clock model on figure 6. In the case q = 5, the transition temperatures T low BKT and T high BKT were treated, like L 0 , as free parameters. The best values were defined as the minimum of the square deviation around the mean over the different strip widths. The collapse on figure 5 was obtained with T low BKT ≃ 0.917 and L 0 ≃ 0.87 at the low temperature BKT transition and T high BKT ≃ 0.927 and L 0 ≃ 0.07 at the second transition. These values are still relatively different with the estimates of T low BKT and T low BKT . However, these values are very sensible to the interval of temperatures on which the mean-square deviation is minimized. Almost as good collapses can be obtained with the temperatures T low BKT ≃ 0.90514(9) and T high BKT ≃ 0.95147(9) [START_REF] Borisenko | Numerical study of the phase transitions in the two-dimensional z(5) vector model[END_REF]. In the case q = 6 plotted on figure 6, we restricted ourselves to use the same values as in Ref. [START_REF] Kumano | Response to a twist in systems with Zp symmetry: The two-dimensional p-state clock model[END_REF].

DMRG study of the Berezinskii-Kosterlitz-Thouless transitions of the 2D five-state clock model9

Fourth-order helicity modulus

As shown by Kim et al. in the case of the XY model [START_REF] Minnhagen | Direct evidence of the discontinuous character of the kosterlitzthouless jump[END_REF], the discontinuity of the helicity modulus Υ 2 at the BKT transition temperature manifests itself as a dip in higher-order derivatives of the free energy density. The temperature dependence of the fourth-order helicity modulus Υ 4 , estimated as [START_REF] Fisher | Helicity modulus, superfluidity, and scaling in isotropic systems[END_REF], is plotted on figure 7. Like the helicity modulus Υ 2 , it vanishes in the high temperature phase. The formation of a dip is clearly seen around T high BKT . However, the low-temperature BKT transition does not manifest itself as a dip. Nevertheless, it seems that the intersection of the curves at temperatures close to T low BKT indicates that Υ 4 takes a universal value, like Υ 2 , at T low BKT rather than displaying a dip. On the right figure, the location of the dip and the temperatures at which the curves intersect for two successive strip widths are plotted versus 1/(ln L) 2 . Both can only be roughly determined but, nevertheless, a linear behavior is observed if one excludes the point L = 16. Linear fits give the BKT transition temperatures T low BKT ≃ 0.91(2) and T high BKT ≃ 0.96(4), again in agreement with earlier estimates T low BKT ≃ 0.90514(9) and T high BKT ≃ 0.95147(9) [START_REF] Borisenko | Numerical study of the phase transitions in the two-dimensional z(5) vector model[END_REF].

Magnetization

Even though the system is infinite in the longitudinal direction of the strip, a finite spontaneous magnetization is observed because magnetic fields are applied at the boundaries. The magnetization m is measured on the site L/2, i.e. at the center of the strip. The data are plotted on figure 8. In the scaling regime L ≫ ξ, m is expected to display an essential singularity ξ -η/2 while, in the finite-size regime, the usual algebraic law with the lattice size, m ∼ L -β/ν with β/ν = η/2, is recovered. The magnetic scaling dimension is usually conjectured to take the same value as in the XY model, i.e. η = 1/4, even though recent Monte Carlo simulations systematically led to smaller values [START_REF] Borisenko | Numerical study of the phase transitions in the two-dimensional z(5) vector model[END_REF]. On figure 8, the scaling hypothesis

m ∼ L -η/2 F m (|T -T BKT | (ln L/L 0 ) 2 ) (11)
is tested with the value η/2 = 1/8. The collapse of the curves corresponding to different strip widths on figure 8 (right) shows that the scaling function F m depends, as expected, only on |T -T BKT | (ln L/L 0 ) 2 . Again, the mean-square deviation was minimized to find to optimal values of the free parameters T low BKT and L 0 . The values T low BKT ≃ 0.917 and L 0 ≃ 0.023 were obtained. As in the case of the helicity modulus, small joint variations of these parameters lead to almost as good collapses so the method is not reliable to determine T low BKT accurately. However, it provides evidence of the scaling behavior [START_REF] Hohenberg | Existence of long-range order in one and two dimensions[END_REF]. We also tried to extract the magnetic exponent η/2 from a fit of magnetization as m ∼ L -η/2 . At the temperature T low BKT ≃ 0.90514, the estimated exponent is 0.116(4), i.e. slightly below the expected value 1/8. However, it increases with temperature, which suggests that strong scaling corrections are present. Using the temperature T low BKT ≃ 0.917 that leads to the best scaling Eq. 11, the exponent takes the value η/2 = 0.1254(30), in perfect agreement with the expected value.

Conclusions

The two BKT phase transitions of the 5-state clock model were studied using the DMRG algorithm. The latter is shown to allow for larger lattice sizes with a better accuracy than Monte Carlo simulations. Despite the Open Boundary Conditions imposed by the technique, the helicity modulus Υ 2 , as well as higher-derivative Υ 4 , can be computed by imposing twisted magnetic fields on the two boundaries of the strip. The topological nature of the two phase transitions is confirmed. Finite-Size displacements of the BKT transition temperatures scale as 1/(ln L) 2 , in agreement with an essential singularity of the correlation length. The exponent σ is therefore the same as for the XY model for both BKT transitions. The helicity modulus is shown to tend towards the universal values predicted by Renormalization Group techniques at the temperatures compatible with the most accurate Monte Carlo estimates. Finally, the scaling behavior of magnetization is also compatible with an essential singularity with the same magnetic scaling dimension η/2 = 1/8 as the XY model.

Figure 1 .

 1 Figure 1. On the left, schematic action of the transfer matrix on a M × L square lattice of the 5-state clock model. The spins on the upper row are denoted as σ 1 to σ L . A new row with spins σ ′ 1 to σ ′ L is added to the system, resulting into a (M + 1) × L lattice. Denote Z M (σ 1 , . . . , σ L ) the partition function of the M × L lattice with fixed spins σ 1 to σ L on the upper row. Similarly, Z M +1 (σ ′ 1 , . . . , σ ′ L ) is the partition function after the addition of the M + 1-th row. The linear relation between them

Figure 2 .

 2 Figure 2. Helicity modulus Υ 2 of the 5-state clock model at the inverse temperature β = 1/k B T = 1.08 versus 1/m 4 where m is the number of states of the truncated vector spaces of the left and right blocks during the DMRG process. The different curves correspond to different strip widths L and the straight lines to linear fits. On the right, a zoom is performed on the curves for L = 128 (top) and L = 256 (bottom).

Figure 3 . 16 T 1 /

 3161 Figure 3. On the left, helicity modulus Υ 2 of the 5-state clock model versus the temperature T for different strip widths L. The two dashed lines correspond to the predicted RG relations at the BKT transition temperatures. On the right, the helicity modulus of the 6-state clock model is plotted for comparison.

Figure 4 .

 4 Figure 4. Temperatures where the helicity modulus Υ 2 of the 5-state clock model intersects the Renormalization Group predictions (8) versus 1/(ln L) 2 . L is the width of the strip. The straight lines are linear fits.

Figure 5 .

 5 Figure 5. Scaling function F Υ (up to a factor) of the 5-state clock model with respect to the scaling argument (β BKT -β) (ln L/L 0 ) 2 . The different symbols correspond to different strip widths. The left figure corresponds to the high-temperature BKT transition while the right one corresponds to the lowtemperature transition.

Figure 6 .

 6 Figure 6. Scaling function F Υ (up to a factor) of the 6-state clock model with respect to the scaling argument (β BKT -β) (ln L/L 0 ) 2 . The different symbols correspond to different strip widths. The left figure corresponds to the high-temperature BKT transition while the right one corresponds to the lowtemperature transition.

16 T 1 /

 161 (ln L)2

Figure 7 .

 7 Figure 7. On the right, fourth derivative Υ 4 of the free energy density (Eq. 7) of the 5-state clock model versus temperature T . The different curves correspond to different strip widths L. The two BKT transitions are signaled by respectively the crossing of the curves and the dip. On the left, temperatures at which the intersection of the curves are observed (red) and location of the dip (black) versus 1/(ln L) 2 . The straight lines are linear fits.

2 mFigure 8 .

 28 Figure 8. On the left, magnetization density m of the 5-state clock model versus the temperature T for different strip widths L. On the right, the scaling function m L 1/8 is plotted with respect to the scaling parameter (β low BKT -β) (ln L/L 0 ) 2 .