Ania Morvannou
email: ania.morvannou@irstea.fr

Nicolas Forquet
email: nicolas.forquet@irstea.fr

Marnik Vanclooster
email: marnik.vanclooster@uclouvain.be

Pascal Molle
email: pascal.molle@irstea.fr

Which Hydraulic Model to Use for Vertical Flow Constructed Wetlands?

Modelling water flow in a Vertical Flow Constructed Wetland (VFWC) is a prerequisite to model wastewater treatment using process based filtering models. The material used in VFCW is very susceptible to generate preferential flow. If this occurs, water will bypass most of the soil porous matrix in a largely unpredictable way. Even if it is possible to simulate water content variations within a VFCW, we cannot correctly model outflow with the standard van Genuchten-Mualem function if preferential flow occurs. A number of various model approaches have been proposed to overcome this problem. These models mostly try to separately describe flow and transport in preferred flow paths and slow or stagnant pore regions. The objective of this study was to study and simulate the hydrodynamic behaviour and the solute transport through a simplified representation of a French VFCW by using both a classical equilibrium model and a non-equilibrium model (dual-porosity model: mobileimmobile water model, with water content mass transfer) included in the HYDRUS-1D software package. Modelling results were compared to a solute breakthrough curve obtained from a tracer experiment, carried out on an existing VFCW. Inlet concentrations were first corrected to take into account tracer loss and experimental uncertainties. Then, both the hydraulic parameters of the mobile and immobile regions ( m r  / m s  and im r  / im s  ) and the transfer coefficient (ω) were optimized to fit the tracer experimental breakthrough curve. The comparison between measured and simulated tracer breakthrough curves indicates that the non-equilibrium approach seem to be the most appropriate for simulating preferential flow paths.

Introduction

Constructed wetlands (CWs) are an attractive wastewater treatment technology for small communities (< 2000 people-equivalents; p.e.), as their simplicity of operation, low cost and reliable treatment efficiency often fit with the limited resources small communities can allocate to wastewater treatment. Irstea (formerly Cemagref) developed an innovative type of first-stage for vertical flow constructed wetlands (VFCW), made of gravel instead of sand. This first-stage type directly accepts raw wastewater, without the need for a preliminary settling tank. Despite good treatment efficiencies, Molle et al. (2008) highlighted that treatment efficiency could be improved if the hydraulic behaviour of the VFCW would be better understood.

The French VFCW are initially entirely composed of gravel but after some years of operation a sludge layer develops on the filter surface and the filter media turns into a matrix made of gravels and organic matter. Water retention is mainly caused by organic matter in the filter. Consequently, both gravitational and capillary flows, due to gravel and the presence of organic matter, respectively, have to be considered (Maier et al., 2009). The macroporosity due to coarse material and the network of reed's stalks may serve as preferential flow paths through which water can bypass most of the porous matrix in a largely unpredictable way. Moreover, during the rest period the sludge dries fast and cracks become visible facilitating preferential flow paths through the sludge layer especially at the beginning of a feeding period (Molle et al., 2006). Consequently, non-equilibrium conditions in pressure heads are created between preferential flow paths and the matrix pore region. Up to now, VFCW is depicted as a porous medium with continuous properties [START_REF] Gerke | A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media[END_REF]Köhne et al., 2009). Variably saturated water flow is then often modelled by the Richards equation assuming that water content and pressure head are at equilibrium at a given node. For simulating preferential flows in structured media, these equilibrium models cannot be applied [START_REF] Gerke | Preferential flow descriptions for structured soils[END_REF]. Moreover, even if it is possible to simulate water content variations within a VFCW using equilibrium models, outflow cannot be correctly modelled with the standard van Genuchten-Mualem function (Morvannou et al, submitted). Therefore, a model representing a non-equilibrium water flow and solute transport would be more appropriate.

Regarding preferential flow models, none has been applied to the flow and transport in VFCW yet. However, lots of studies have been carried out concerning preferential flows in the field of hydrology for simulating water flow and/or solute transport in saturated and unsaturated porous soils through macropores (Šimůnek et al., 2001;Haws et al., 2005;Köhne et al., 2006;Maier et al., 2009). Šimunek et al., (2003, 2008) carried out extensive literature reviews of existing models. Most of them require the parallel consideration of two or more domains for the same porous medium. Four of them have been implemented in the HYDRUS software package (Šimůnek et al., 2006). These models have in common that water flow is simulated using the Richards equation and solute transport using the convection-dispersion equation. In contrast, they differ by the way flow and transport in ''fast" (macropore domain) and slow or stagnant pore space (matrix domain) are modelled. In addition they are coupled by kinetic mass transfer terms between the two pore domains. However, the disadvantage of these approaches is that the computational demand significantly increases as well as the number of additional parameters, which are hard to estimate (Šimůnek and van Genuchten, 2008).

The purpose of this paper is to study and simulate the hydrodynamic behaviour and the solute transport through a simplified representation of a VFCW by using both a classical equilibrium model and a non-equilibrium model included in the HYDRUS-1D software package. We used the dual-porosity model as the simplest non-equilibrium approach, with fewer parameters, while having both non-equilibrium flow and solute transport. Modelling results will be compared to a solute breakthrough curve obtained from a tracer experiment carried out on an existing VFCW to determine which approach is the most appropriate. Finally, using a more realistic approach we tried to calibrate the dual-porosity model.

Materials and methods

Tracer experiment

The Evieu wastewater treatment plant (200 p.e., Ain, France) is in operation since 2004. The first stage is made up of three VFCWs receiving raw wastewater. The pump sump at the plant inlet only performs mechanical screening using a 5 cm mesh screen. Each filter is fed according to a feeding/rest regime of 3.5/7 days. Effluent from this first stage is then connected to a second pump sump and separated between the second stage's vertical and horizontal flow constructed wetlands according to the experimental objectives (Molle et al., 2008). All beds are planted with Phragmites australis.

A fluorescein infiltration experiment was carried out from 04/15/2010 to 04/19/2010, corresponding to the entire feeding period, on one of the VFCW of the first stage of the wastewater treatment plant. The selected vertical filter (28 m², 2.9 m wide x 9.7 m long) is designed according to French recommendations. From the bottom to the top, it contains a 15 cm thick drainage layer (grain size of 30-60 mm), a 10 cm transition layer (grain size 15-25 mm) and a 60 cm gravel layer (d 10 = 2.46 mm; UC = d 60 /d 10 = 1.39; average porosity of 40.4%). As the filter has been working at nominal load for 7 years, a sludge layer of about 20 cm has developed at the top of the filter. The VFCW is fed by raw wastewater in batches of 5 cm at a rate of 1.23 m 3 /h/m² on average. During spells of dry weather, 8 separate batches are processed per day. Water is drained by a 160 mm diameter drainage pipe (0.42 m of pipe per m²), allowing passive aeration from the bottom as well.

For the tracer experiment, 1.98 g of fluorescein was diluted in a batch volume (1973 L) and applied to the filter. This is equal to 0.07 g/m² and the applied solution fluorescein concentration was 0.001 g/L. It should be noted that all the fluorescein was not applied in only one batch as expected; some remained in the pump sump (despite we emptied it after the tracer batch) and was consequently applied with the next batches. Unfortunately, the remaining concentration could not be determined. In addition, after the batch some fluorescein remains also on the top of the filter due to the low permeability and the heterogeneity of the sludge deposit. After the fluorescein application, its concentration was recorded at the VFCW outlet (online fluorimeter PELI 1300 Case, Peli Products, USA) along with the outflow rate (measured by a venture flume) at a time interval of 1 minute throughout the tracer experiment.

At the same time, water contents were recorded with TDR probes, but only during the first batch of the tracer experiment. TDR probes were inserted at three different depths of the filter (13, 31 and 47 cm below the VFCW top surface) and connected through multiplexers to a signal generator/analyzer (model TDR100, Campbell ScientificTM, Logan, UT). Results were recorded using a datalogger (model CR1000, Campbell ScientificTM Logan, UT). Measurements were made at 20 second time intervals.

An introduction to dual-porosity models

The aim of this section is not to provide a full mathematical description of the dual-porosity models. For this purpose the reader may refer to Šimůnek et al. (2003, 2008). Our objectives are to emphasize how dual-porosity models represent preferential flows and to highlight how complexity is increased by additional parameters. According to Šimůnek et al. (2003), preferential flows have two main characteristics: (1) the ability to quickly propagate by "bypassing a large part of the matrix pore-space", and (2) the non-equilibrium between the water pressure head within the macropores and the water content in the rest of the matrix. Non-equilibrium is a key concept of preferential flow modelling. In the classical description of flow in variably saturated porous media based on the Richards equation, it is stated that the water pressure head and water content are always at equilibrium and linked by a water retention relationship (e.g. the van Genuchten relationship). This assumption does not hold anymore in case preferential flow occurs.

The dual-porosity model does not aim at physically representing the paths of quick flow (it would be very challenging to map every macropore at the scale of a VFCW) but rather to incorporate their effects in the representative elementary volume (REV) that serves as a basis for modelling. Figure 1 comparatively represents how macropores can be observed on the filter and how they are modelled at the REV scale. (1) In dual-porosity models, it is assumed that flow only takes place within the macropores (the mobile region), whereas the flow in the rest of the matrix is neglected. Water in the rest of the matrix constitutes the immobile region. The hydrodynamic parameters (α, n and K s ) are representative of the flow in macropores and therefore must be adapted, especially K s , to produce quick outflow. The total porosity is divided between the mobile ( m s  ) and the immobile ( im s  ) regions. This brings additional complexity to the model: for a given saturated hydraulic conductivity, the larger is the mobile region, the slower is the flow. For sake of simplicity we considered residual water content to be null for both regions.

(2) Non-equilibrium is modelled by adding a time-dependent parameter (ω) to the relationship between water content in the immobile region and water content in the mobile region that is linked to the water pressure head. This is a very important feature that allows water exchange at different time scale than the main gravitational drainage. This parameter also represents an increase in model complexity and needs to be calibrated for each material in the model.

Model description

Based on the significant ponding occuring during feeding (~5 cm/m²), we assumed that water distribution is homogeneous. Consequently, we modelled the water flow and the solute transport in one dimension with HYDRUS-1D. Moreover, in contrast to the previous study carried out on the same VFCW where we considered four layers (the sludge layer, the first gravel-sludge layer, the second gravel-sludge layer, and a final practically gravel-only layer), in this study we restrained the domain to only two layers: the sludge and the gravel-only layer. As the number of parameters in the dual-porosity model is higher than in the equilibrium model, a simplified representation of the VFCW facilitates modelling and understanding of the non-equilibrium phenomenon.

The one-dimensional mesh used for simulations consisted of 101 equidistant nodes. The top boundary condition is a time-dependent atmospheric boundary condition with a threshold value for surface runoff set to 10 cm. Wastewater load duration and flow rate are specified, while evapotranspiration is neglected. If the incoming flow rate exceeds infiltration capacity, the ponding of water above the surface is then taken into account until it reaches 10 cm, a value that was never observed in our case. The bottom boundary condition is a seepage face. Initial conditions are set to hydrostatic equilibrium.

Tested hydraulic models

Water flow and solute transport were simulated using two different modelling approaches: the generally-equilibrium model and the dual-porosity model (non-equilibrium model). Two dual-porosity models are available in HYDRUS-1D depending on if the mass transfer rate (Γ w ) is proportional to the difference in effective water contents, or to the difference in effective pressure heads [START_REF] Gerke | A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media[END_REF]. We used the first model because it requires significantly fewer parameters since one does not need to know the retention function for the matrix region explicitly, but only its residual and saturated water contents.

                    im r im s im r im m r m s m r m im e m e im w S S t            (1)
Which, in our case, simplifies to:

         im s im m s m w      ( 2 )
where Γ w is the mass transfer rate for water between the mobile and immobile regions, θ m and θ im are the mobile and immobile water content, respectively, m Even if the second model provides a more realistic description of the exchange rate between the fracture and matrix regions (the difference in pressure heads is the actual driving force for mass transfer), the mass transfer term may be inherently more unstable numerically (Šimůnek et al., 2003). For the equilibrium model, we used hydraulic parameters determined from a previous study (Morvannou, 2012) 1).

The first aim of this study was to have a better insight into hydrodynamic and solute transport phenomena taking place inside a VFCW. In a first attempt we simulated the tracer application with only one pulse. Moreover, in order to evaluate if the observed multiple tracer peaks could be due to tracer kept in the sump or at the surface of the filter, we also tested a case with seven pulses of decreasing concentration. Modelling results from the equilibrium and non-equilibrium models were compared in terms of concentrations and cumulative fluxes computed at the outlet.

Calibration of the dual-porosity model

The second objective was to calibrate the dual-porosity model. In this part of the study, we divided the gravel-sludge layer in two layers in order to take into account the decreasing gradient of organic matter content versus depth and to compare the results with the water content measured by the TDR probes installed in the filter. The one-dimensional mesh, the top and bottom boundary conditions were the same as those used for the previous simulations. The setup of initial conditions necessitated a prior initialization step. We could not use measured water contents (TDR probes) since they do not provide continuous pressure head profile at the initial time and cause the model not to converge (Radcliffe and Šimůnek, 2010). The prior initialization step consists in repeating a batch simulation until pressure heads reach a pseudo-permanent state. We then selected the time for which water content profile matches the experimentally-measured initial water content. The corresponding water pressure profile was then considered as initial condition. Inlet tracer concentrations were adapted from the tracer concentrations measured at the outlet. Therefore, we applied four batches with decreasing tracer concentrations. Then, both the hydraulic parameters of the mobile and immobile regions ( m r  / m s  and im r  / im s  ) and the transfer coefficient (ω) were optimized to fit the tracer experimental breakthrough curve. The total tracer recovery is 62% and consequently reveals some fluorescein losses. They could be due to the ponding of water above the surface (photodegradation of the tracer), tracer retention on the organic matter (adsorption) and uncertainties of flow and concentration measurements.

Results and discussion

Tracer experiment

Figure 2 highlights two main phenomena. (1) Not only one but several peaks are observed on the experimental breakthrough curve (Figure 2 (1)). Even if the tracer is mixed only once with the incoming wastewater, it seems like the tracer application was carried out in several times (remaining tracer in the sump and on filter surface). The tracer leaves the filter from the first batch pointing out the presence of preferential flows in the filter. (2) The second phenomenon is a change of the breakthrough curve shape. After the seven first peaks, the tracer concentration decreases when a new batch arrives and then increases according with a logarithmic tendency towards an equilibrium value (Figure 2 (2)).

Modelling results

Table 1 presents the hydrodynamic parameter values used for the equilibrium and nonequilibrium models.

Table 1. Parameter values of the hydrodynamic properties of the column representing the VFCW used in HYDRUS-1D When comparing the simulated tracer breakthrough curve obtained from the two models, we observe that the tracer leaves the filter faster with the dual-porosity model than with the equilibrium one. On the other end, we note that the equilibrium model is not able to reproduce the shape inversion observed with the experimental tracer breakthrough curve (Figure 2 (2)). On the contrary, the dual-porosity provides the same shape inversion due to the tracer transfer into the immobile region and its release after each batch. Preferential flow effects modelled as mobile and immobile water are very important for matching fluorescein tracer transport characteristics.

Equilibrium model

Layer

θ r [-] θ s [-] α [1/cm] n [-] K s [cm/s] λ [-] 1 0.
 [-] m s  [-] α m [1/cm] n m [-] K s [cm/s] λ [-] im r  [-] im s  [-] ω [1/s]
Secondly, applying only one batch containing tracer does not allow simulating the multiple peaks of tracer observed in Figure 2. It is necessary to apply the tracer several times to reproduce the peaks of tracer observed at the beginning of the tracer experiment.

If we just compare the cumulative tracer fluxes in the outlet of the filter we cannot see an important difference between the equilibrium and non-equilibrium results. Indeed, cumulative fluxes are not sensitive enough for comparing the modelling results whereas tracer concentrations provide information about non-equilibrium phenomenon which takes place in the filter.

The comparison between measured and simulated tracer breakthrough curves indicates that the non-equilibrium approach seem to be the most appropriate for simulating preferential flow paths.

Calibration results

Table 2 present the hydraulic parameter values for the dual-porosity model adapted from both the tracer concentrations obtained at the outlet of the filter and the water contents recorded in the VFCW. After adapting hydraulic parameters, the simulation using the dual-porosity model (Figure 4) showed a good match with the measured data (RMSE = 1.47.10 -5 ). We identified the four pulses corresponding to the four inlet concentrations at the beginning of the tracer experiment. After, these four batches, we observed the same tendency that we observed in the Figure 2 (2): the tracer concentration decreases right after the loading and then increases according with a logarithmic tendency towards an equilibrium value until a new batch.

m r  [-] m s  [-] α m [1/cm] n m [-] K s [cm/s] λ [-] im r  [-] im s  [-] ω [1/
First, we saw that the water retention and hence the water content decreases with depth, as organic matter content is lower deeper in the filter. Secondly, even if we only recorded water contents for the first tracer batch, Figure 5 showed that the simulation of water contents was good for the sludge layer (13 cm) where predicted water contents were close to observed water contents (RMSE = 9.1.10 -3 ). The simulated water contents for the first (31 cm) and second (47 cm) gravel-sludge layers were moderately good (RMSE = 0.39, 0.43, respectively).

Therefore, the dual-porosity model is able to simulate both the tracer breakthrough and water contents at different depths of the VFCW.

Conclusions

A comparison was carried out between a classical equilibrium model and a non-equilibrium model (dual-porosity model) for simulating the hydrodynamic behaviour and the solute transport through a simplified representation of a French VFCW. The HYDRUS-1D software package was used for modelling flow and transport. Modelling results were compared to a solute breakthrough curve obtained from a tracer experiment carried out on an existing VFCW. This allows the determination of the most appropriate modelling approach.

The comparison between measured and simulated tracer breakthrough curves indicated that the tracer leaves the filter faster with the dual-porosity model than with the equilibrium one, which is consistent with the experimental tracer breakthrough. Moreover, the dual-porosity model provides the same shape inversion due to the tracer transfer into the immobile region and its release after each batch. Moreover, in addition to simulate tracer breakthrough, the dual-porosity model succeeds to simulate the water content profiles at different depths of the filter. Therefore, the non-equilibrium approach seems to be the most appropriate for simulating preferential flow paths in a French-type VFCW.

Figure 1 .

 1 Figure 1. Schematic representation of macropores and preferential flow paths in a VFCW (left) and the corresponding conceptual physical non-equilibrium model for water flow and solute transport (right)



  are the mobile and immobile saturated water content, respectively, ω is a first-order rate coefficient (T -1 ), and m e S and im e S are effective fluid saturations of the mobile and immobile regions, respectively.

  and for the dual-porosity model some hydraulic parameter values were adapted from Morvannou (2012) ( m s  , α m , n m , λ and im s  ) whereas others were arbitrarily fixed ( m r  , K s , im r  and ω) (Table

Figure 2

 2 Figure 2 presents the tracer concentration measured at the outlet of the VFCW and the cumulative percentage of the tracer recovery.

Figure 2 .

 2 Figure 2. Tracer experimental breakthrough curve and cumulative tracer recovery percentage

Figure 3

 3 Figure 3 presents the simulated concentrations and fluxes obtained from the equilibrium and non-equilibrium models. It also shows the results obtained with only one tracer pulse and those with seven tracer pulses.

Figure 3 .

 3 Figure 3. Simulated concentrations (left) and cumulative fluxes (right) for one pulse (top) or seven pulses (bottom) from the equilibrium and non-equilibrium models

Figure 4 .

 4 Figure 4. Measured (black line) and simulated (red line) concentrations from the tracer experiment and the dual-porosity model, respectively

Table 2 .

 2 Parameter values used for the calibration of the dual-porosity model representing the VFCW behaviour in HYDRUS-1D

	Layer

4th International Conference HYDRUS Software Applications to SubsurfaceFlow and Contaminant Transport Problems, Prague,[21][22]

Haws, N.W., P.S.C. Rao, J. Šimůnek, and I.C. Poyer, Single-porosity and dual-porosity modeling of water flow and solute transport in subsurface-drained fields using effective field-scale parameters, Journal of Hydrology, 313(3-4), 257-273, 2005. Köhne, J.M., S. Köhne, and J. Šimůnek, Multi-process herbicide transport in structured soil columns: experiments and model analysis, Journal of Contaminant Hydrology, 85, 1-32, 2006. Köhne, J.M., S. Köhne, and J. Šimůnek, A review of model applications for structured soils: a) Water flow and tracer transport, Journal of Contaminant Hydrology, 104, 4-35, 2009. Maier, U., C. DeBiase, O. Baeder-Bederski, and P. Bayer, Calibration of hydraulic parameters for large-scale vertical flow constructed wetlands, Journal of Hydrology, 369(3-4), 260-273, 2009. Molle, P., A. Liénard, A. Grasmick, and A. Iwema, Effect of reeds and feeding operations on hydraulic behavior of vertical flow constructed wetlands under hydraulic overloads, Water Research, 40(3), 606-612, 2006. Molle, P., S. Prost-Boucle, and A. Liénard, Potential for total nitrogen removal by combining vertical flow and horizontal flow constructed wetlands: a full-scale experiment study, Ecological Engineering, 34(1), 23-29, 2008. Morvannou, A., Dynamic modelling of nitrification in vertical flow constructed wetlands. Sciences agronomiques et ingénierie biologique, Université catholique de Louvain, 184p, 2012. Morvannou, A., N. Forquet, M. Vanclooster, and P. Molle (Submitted), On the characterization of the hydraulic properties of the filter material of a Vertical Flow Constructed Wetland. Radcliffe, D., and J. Šimůnek, Soil Physics with HYDRUS: Modelling and Applications, CRC Press, Taylor and Francis Group, 373 p, 2010. Šimůnek, J., O. Wendroth, N. Wypler, and M. Th. van Genuchten, Non-equilibrium water flow characterized by means of upward infiltration experiments, European Journal of Soil Science, 52, [13][14][15][16][17][18][19][20][21][22][23][24] 2001. Šimůnek, J., N.J. Jarvis, M. Th. van Genuchten, and A. Gärdenäs, Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone, Journal of Hydrology, 272(1-4), 2003. Šimůnek, J., D. Jacques, M. Th. van Genuchten, and D. Mallants, Multicomponent geochemical transport modelling using the HYDRUS computer software packages, Journal America Water Resources Association, 42(6), 1537-1547, 2006. Šimůnek, J., and van Genuchten, M. Th. (2008). Modeling nonequilibrium flow and transport with HYDRUS, Vadose Zone Journal, 7(2), 782-797 Special Issue "Vadose Zone Modeling", 2008.


Jean-Luc Fouquet

Henri Thuillier

Jean-Marie Vanherpe

On a family of cubic graphs containing the ower snarks

Keywords: cubic graph, perfect matching, strong matching, counting, hamiltonian cycle, 2-factor hamiltonian

We consider cubic graphs formed with k ≥ 2 disjoint claws C i ∼ K 1,3 (0 ≤ i ≤ k -1) such that for every integer i modulo k the three vertices of degree 1 of C i are joined to the three vertices of degree 1 of C i-1 and joined to the three vertices of degree 1 of C i+1 . Denote by t i the vertex of degree 3 of C i and by T the set {t 1 , t 2 , ..., t k-1 }.

). A strong matching M in a graph G is a matching M such that there is no edge of E(G) connecting any two edges of M . A cubic graph having a perfect matching union of two strong matchings is said to be a Jaeger's graph. We characterize the graphs F S(j, k) that are Jaeger's graphs.

Introduction

The complete bipartite graph K 1,3 is called, as usually, a claw. Let k be an integer ≥ 2 and let G be a cubic graph on 4k vertices formed with k disjoint claws C i = {x i , y i , z i , t i } (0 ≤ i ≤ k -1) where t i (the center of C i ) is joined to the three independent vertices x i , y i and z i (the external vertices of C i ). For every integer i modulo k C i has three neighbours in C i-1 and three neighbours in C i+1 . For any integer k ≥ 2 we shall denote the set of integers modulo k as Z k . In the sequel of this paper indices i of claws C i belong to Z k . By renaming some external vertices of claws we can suppose, without loss of generality, that {x i x i+1 , y i y i+1 , z i z i+1 } are edges for any i distinct from k -1. That is to say the subgraph induced on X = {x 0 , x 1 , . . . , x k-1 } (respectively Y = {y 0 , y 1 , . . . , y k-1 }, Z = {z 0 , z 1 , . . . , z k-1 }) is a path or a cycle (as induced subgraph of G). Denote by T the set of the internal vertices {t 0 , t 1 , . . . , t k-1 }.

Up to isomorphism, the matching joining the external vertices of C k-1 to those of C 0 (also called, for k ≥ 3, edges between C k-1 and C 0 ) determines the graph G. In this way we construct essentially three distinct graphs, namely F S(1, k), F S(2, k) and F S(3, k). The graph F S(j, k) (j ∈ {1, 2, 3}) is the graph where the set of vertices ∪ i=k-1 i=0 {C i \ {t i }} induces j cycles. For k ≥ 3 and any j ∈ {1, 2, 3} the graph F S(j, k) is a simple cubic graph. When k is odd, the F S (2, k) are the graphs known as the ower snarks [START_REF] Isaacs | Innite families of non-trivial trivalent graphs which are not Tait colorable[END_REF]. We note that F S(3, 2) and F S(2, 2) are multigraphs, and that F S(1, 2) is isomorphic to the cube. For k = 2 the notion of "edge between C k-1 and C 0 " is ambiguous, so we must dene it precisely. For two parallel edges having one end in C 0 and the other in C 1 , for instance two parallel edges having x 0 and x 1 as endvertices, we denote one edge by x 0 x 1 and the other by x 1 x 0 . An edge in {x 1 x 0 , x 1 y 0 , x 1 z 0 , y 1 x 0 , y 1 y 0 , y 1 z 0 , z 1 x 0 , z 1 y 0 , z 1 z 0 }, if it exists, is an edge between C 1 and C 0 . We will say that x 0 x 1 , y 0 y 1 and z 0 z 1 are edges between C 0 and C 1 .

By using an ad hoc translation of the indices of claws (and of their vertices) and renaming some external vertices of claws, we see that for any reasoning about a sequence of h ≥ 3 consecutive claws (C i , C i+1 , C i+2 , . . . , C i+h-1 ) there is no loss of generality to suppose that 0 ≤ i < i+h-1 ≤ k -1. For a sequence of claws (C p , ..., C r ) with 0 ≤ p < r ≤ k -1, since 0 is a possible value for subscript p and since k -1 is a possible value for subscript r, it will be useful from time to time to denote by x ′ p-1 the neighbour in C p-1 of the vertex x p of C p (recall that x ′ p-1 ∈ {x k-1 , y k-1 , z k-1 } if p = 0), and to denote by x ′ r+1 the neighbour in C r+1 of the vertex x r of C r (recall that x ′ r+1 ∈ {x 0 , y 0 , z 0 } if r = k -1). We shall make use of analogous notations for neighbours of y p , z p , y r and z r .

We shall prove in the following lemma that there are essentially two types of perfect matchings in F S(j, k).

Lemma 1 Let G ∈ {F S(j, k), j ∈ {1, 2, 3}, k ≥ 2} and let M be a perfect matching of G. Then the 2-factor G \ M induces a path of length 2 and an isolated vertex in each claw C i (i ∈ Z k ) and M fulls one (and only one) of the three following properties : i) For every i in Z k M contains exactly one edge joining the claw C i to the claw C i+1 , ii) For every even i in Z k M contains exactly two edges between C i and C i+1 and none between C i-1 and C i , iii) For every odd i in Z k M contains exactly two edges between C i and C i+1 and none between C i-1 and C i .

Moreover, when k is odd M satises only item i).

Proof Let M be a perfect matching of G = F S(j, k) for some j ∈ {1, 2, 3}.

Since M contains exactly one edge of each claw, it is obvious that G \ M induces a path of length 2 and an isolated vertex in each claw C i .

For each claw C i of G the vertex t i must be saturated by an edge of M whose end (distinct from t i ) is in {x i , y i , z i }. Hence there are exactly two edges of M having one end in C i and the other in

C i-1 ∪ C i+1 .
If there are two edges of M between C i and C i+1 then there is no edge of M between C i-1 and C i . If there are two edges of M between C i-1 and C i then there is no edge of M between C i and C i+1 . Hence, we get ii) or iii) and we must have an even number k of claws in G.

Assume now that there is only one edge of M between C i-1 and C i . Then there exists exactly one edge between C i and C i+1 and, extending this trick to each claw of G, we get i) when k is even or odd.

Denition 2 We say that a perfect matching M of F S(j, k) is of type 1 in Case i) of Lemma 1 and of type 2 in Cases ii) and iii). If neccessary, to distinguish Case ii) from Case iii) we shall say type 2.0 in Case ii) and type 2.1

in Case iii). We note that the numbers of perfect matchings of type 2.0 and of type 2.1 are equal.

Notation : The length of a path P (respectively a cycle Γ) is denoted by l(P )

(respectively l(Γ)).

2 Counting perfect matchings of F S(j, k)

We shall say that a vertex v of a cubic graph G is inated into a triangle when we construct a new cubic graph G ′ by deleting v and adding three new vertices inducing a triangle and joining each vertex of the neighbourhood N (v) of v to a single vertex of this new triangle. We say also that G ′ is obtained from G by a triangular extension. The converse operation is the contraction or reduction of the triangle. The number of perfect matchings of G is denoted by µ(G).

Lemma 3 Let G be a bipartite cubic graph and let {V 1 , V 2 } be the bipartition of its vertex set. Assume that each vertex in some subset W 1 ⊆ V 1 is inated into a triangle and let G ′ be the graph obtained in that way. Then

µ(G) = µ(G ′ ).
Proof Note that {V 1 , V 2 } is a balanced bipartition and, by König's Theorem, the graph G is a cubic 3-edge colourable graph. So, G ′ is also a cubic 3-edge colourable graph (hence, G and G ′ have perfect matchings). Let M be

a perfect matching of G ′ . Each vertex of V 1 \ W 1 is saturated by an edge whose second end vertex is in V 2 . Let A ⊆ V 2 be the set of vertices so saturated in V 2 .
Assume that some triangle of G ′ is such that the three vertices are saturated by three edges having one end in the triangle and the second one in V 2 . Then we need to have at least |W 1 | + 2 vertices in V 2 \ A, a contradiction. Hence, M must have exactly one edge in each triangle and the contraction of each triangle in order to get back G transforms M in a perfect matching of G.

Conversely, each perfect matching of G leads to a unique perfect matching of G ′ and we obtain the result.

Let us denote by µ(j, k) the number of perfect matchings of F S(j, k), µ 1 (j, k) its number of perfect matchings of type 1 and µ 2 (j, k) its number of perfect matchings of type 2.

Lemma 4 We have

• µ(1, 3) = µ 1 (1, 3) = 9 • µ(2, 3) = µ 1 (2, 3) = 8 • µ(3, 3) = µ 1 (3, 3) = 6 • µ(1, 2) = 9, µ 1 (1, 2) = 3 • µ(2, 2) = 10, µ 1 (2, 2) = 4 • µ(3, 2) = 12, µ 1 (3, 2) = 6
Proof The cycle containing the external vertices of the claws of the graph

F S(1, 3) is x 0 , x 1 , x 2 , y 0 , y 1 , y 2 , z 0 , z 1 , z 2 , x 0 . Consider a perfect matching M
containing the edge t 0 x 0 . There are two cases: i) x 1 x 2 ∈ M and ii) x 1 t 1 ∈ M . In Case i) we must have y 0 y 1 , t 1 z 1 , t 2 z 2 , z 0 y 2 ∈ M . In Case ii) there are two sub-cases: ii).a x 2 y 0 ∈ M and ii).b x 2 t 2 ∈ M . In Case ii).a we must have y 1 y 2 , t 2 z 2 , z 0 z 1 ∈ M and in Case ii).b we must have y 0 y 1 , y 2 z 0 , z 1 z 2 ∈ M . Thus, there are exactly 3 distinct perfect matching containing t 0 x 0 . By symmetry, there are 3 distinct perfect matchings containing t 0 y 0 , and 3 distinct matchings containing t 0 z 0 , therefore µ(1, 3) = 9.

It is well known that the Petersen graph has exactly 6 perfect matchings. Since F S(2, 3) is obtained from the Petersen graph by inating a vertex into a triangle these 6 perfect matchings lead to 6 perfect matchings of F S(2, 3). We have two new perfect matchings when considering the three edges connected to this triangle (we have two ways to include these edges into a perfect matching). Hence µ(2, 3) = 8.

F S(3, 3) is obtained from K 3,3 by inating three vertices in the same colour of the bipartition. Since K 3,3 has six perfect matchings, applying Lemma 3 we get immediately the result for µ(3, 3).

Is is a routine matter to obtain the values for F S(j, 2) (j ∈ {1, 2, 3}).

Theorem 5 The numbers µ 1 (j, k) of perfect matchings of type 1 of F S(j, k)

(j ∈ {1, 2, 3}
) are given by:

• µ 1 (1, k) = 2 k -(-1) k • µ 1 (2, k) = 2 k • µ 1 (3, k) = 2 k + 2(-1) k
Proof By Lemma 4, when k = 2 or k = 3 µ 1 (j, k) fulls our requirement. We want to compute µ 1 (j, k) by induction on k and we assume that the property holds for F S(j, k -2) with k ≥ 4 (j ∈ {1, 2, 3}).

The following trick will be helpful. Let us consider the four consecutive claws C k-2 , C k-1 , C 0 and C 1 of F S(j, k) (j ∈ {1, 2, 3}). We can delete C k-1 and C 0 and join the three external vertices of C k-2 to the three external vertices of C 1 by a matching (to choose) in such a way that the resulting graph is a graph F S(j ′ , k -2), with j ′ ∈ {1, 2, 3} (it must be clear that j and j ′ may be dierent). In order to count the number of perfect matchings of type 1 of F S(j, k) we need to consider nine numbers, each of them counts the number of perfect matchings of type 1 of F S(j, k) using two edges, one being an edge between C k-2 and C k-1 and the other being an edge between C 1 and C 0 .

Let us dene ν(e, e ′ ) as the number of perfect matchings of type 1 containing the two edges e and e ′ . Then we set

a 1 = ν(x k-2 x k-1 , x 0 x 1 ), a 2 = ν(x k-2 x k-1 , y 0 y 1 ), a 3 = ν(x k-2 x k-1 , z 0 z 1 ) a 4 = ν(y k-2 y k-1 , x 0 x 1 ), a 5 = ν(y k-2 y k-1 , y 0 y 1 ), a 6 = ν(y k-2 y k-1 , z 0 z 1 ) a 7 = ν(z k-2 z k-1 , x 0 x 1 ), a 8 = ν(z k-2 z k-1 , y 0 y 1 ), a 9 = ν(z k-2 z k-1 , z 0 z 1 )
Obviously we have µ 1 (j, k) = 9 i=1 a i . Let us delete the vertices of C k-1 and C 0 and denote by H the resulting induced subgraph.

Claim 1 When j = 1, we have

a 2 + a 6 + a 7 = 2µ 1 (1, k -2) a 1 + a 5 + a 9 = µ 1 (3, k -2) a 3 + a 4 + a 8 = µ 1 (1, k -2)
Proof Without loss of generality we can consider that x k-1 y 0 , y k-1 z 0 and z k-1 x 0 are edges of F S(1, k).

In order to evaluate a 2 + a 6 + a 7 we add the edges x k-2 y 1 , y k-2 z 1 and z k-2 x 1 to H. In other words we set

H 1 = (V (H), E(H) ∪ {x k-2 y 1 , y k-2 z 1 , z k-2 x 1 }), we get hence F S(1, k -2)
. Each perfect matching of type 1 of H 1 containing x k-2 y 1 gives two perfect matchings of type 1 of F S(1, k) as well as each perfect matching of type 1 of H 1 containing y k-2 z 1 and each perfect matching of type

1 containing the edge z k-2 x 1 . If follows that a 2 + a 6 + a 7 = 2µ 1 (1, k -2). Let us now set H 2 = (V (H), E(H)∪{x k-2 x 1 , y k-2 y 1 , z k-2 z 1 }). The graph H 2 is isomorphic to F S(3, k -2). Each perfect matching of type 1 of H 2 containing x k-2 x 1 gives one perfect matching of type 1 of F S(1, k) as well as each perfect matching of type 1 of H 2 containing y k-2 y 1 and each perfect matching of type 1 of H 2 containing z k-2 z 1 . Consequently a 1 + a 5 + a 9 = µ 1 (3, k -2).
For computing a 3 +a 4 +a 8 we set

H 3 = (V (H), E(H)∪{x k-2 z 1 , y k-2 x 1 , z k-2 y 1 }), that is a graph isomorphic to F S(1, k -2). Each perfect matching of type 1 of H 3 containing x k-2 z 1 gives one perfect matching of type 1 of F S(1, k) as well
as each perfect matching of type 1 of H 3 containing y k-2 x 1 and each perfect matching of type

1 of H 3 containing z k-2 y 1 . Thus a 3 + a 4 + a 8 = µ 1 (1, k -2) as claimed.
Claim 2 When j = 2, we have

a 1 + a 6 + a 8 = 2µ 1 (2, k -2) a 2 + a 4 + a 9 = µ 1 (2, k -2) a 3 + a 5 + a 7 = µ 1 (2, k -2)
Proof Without loss of generality we can consider that x k-1 x 0 , y k-1 z 0 and z k-1 y 0 are edges of F S(2, k).

In order to evaluate a 1 + a 6 + a 8 we add the edges

x k-2 x 1 , y k-2 z 1 and z k-2 y 1 to H. In other words we set H 1 = (V (H), E(H) ∪ {x k-2 x 1 , y k-2 z 1 , z k-2 y 1 }), we get hence F S(2, k-2). Each perfect matching of type 1 of H 1 leads to precisely 2 perfect matchings of type 1 of F S(2, k), thus a 1 + a 6 + a 8 = 2µ 1 (2, k -2).
We consider now the graph

H 2 = (V (H), E(H) ∪ {x k-2 y 1 , y k-2 x 1 , z k-2 z 1 }) isomorphic to F S(2, k -2). Each perfect matching of type 1 of H 2 gives one perfect matching of type 1 of F S(2, k), consequently a 2 +a 4 +a 9 = µ 1 (2, k -2). Let H 3 = (V (H), E(H) ∪ {x k-2 z 1 , y k-2 y 1 , z k-2 x 1 }), a graph isomorphic to F S(2, k -2). Each perfect matching of type 1 of H 3 leads to precisely one per- fect matching of type 1 of F S(2, k), consequently a 3 + a 5 + a 7 = µ 1 (2, k -2).
Claim 3 When j = 3, we have

a 1 + a 5 + a 9 = 2µ 1 (3, k -2) a 2 + a 6 + a 7 = µ 1 (1, k -2) a 3 + a 4 + a 8 = µ 1 (1, k -2) Proof Suppose that x k-1 x 0 , y k-1 y 0 , z k-1 z 0 are edges of F S(3, k).
In order to evaluate a 1 + a 5 + a 9 we add the edges x k-2 x 1 , y k-2 y 1 and z k-2 z 1 to H. In other words we set

H 1 = (V (H), E(H) ∪ {x k-2 x 1 , y k-2 y 1 , z k-2 z 1 }), we get hence F S(3, k -2). Each perfect matching of type 1 of H 1 can be extended into 2 perfect matchings of type 1 of F S(3, k), thus a 1 +a 5 +a 9 = 2µ 1 (3, k -2). Let H 2 be isomorphic to F S(1, k -2) with the edges x k-2 y 1 , y k-2 z 1 , z k-2 x 1 . A perfect matching of type 1 of H 2 can precisely be extended into one perfect matching of type 1 of F S(3, k), then a 2 + a 6 + a 7 = µ 1 (1, k -2). Let H 3 = (V (H), E(H) ∪ {x k-2 z 1 , y k-2 x 1 , z k-2 y 1 }), a graph isomorphic to F S(1, k -2). Each perfect matching of type 1 of H 3 gives one perfect match- ing of type 1 of F S(3, k), consequently a 3 + a 4 + a 8 = µ 1 (1, k -2).
Recall that µ 1 (j, k) = 9 i=1 a i . Then it follows from Claims 1, 2 and 3 :

µ 1 (1, k) = 3µ 1 (1, k -2) + µ 1 (3, k -2) µ 1 (2, k) = 4µ 1 (2, k -2) µ 1 (3, k) = 2µ 1 (3, k -2) + 2µ 1 (1, k -2)
By induction,

µ 1 (1, k -2) = 2 k-2 -(-1) k µ 1 (2, k -2) = 2 k-2 µ 1 (3, k -2) = 2 k-2 + 2(-1) k Thus , µ 1 (1, k) = 3(2 k-2 -(-1) k ) + 2 k-2 + 2(-1) k = 2 k -(-1) k µ 1 (2, k) = 4(2 k-2 ) = 2 k µ 1 (3, k) = 2(2 k-2 + 2(-1) k ) + 2(2 k-2 -(-1) k ) = 2 k + 2(-1) k
This ends the proof.

Theorem 6 The number µ 2 (j, k) of perfect matchings of F S(j, k) (j ∈ {1, 2, 3}) veries :

µ 2 (j, k) = 0 when k is odd and µ 2 (j, k) = 2 × 3 k 2 otherwise.
Proof When k is odd, we have µ 2 (j, k) = 0 by Lemma 1.

When k is even, let M be a perfect matching of type 2 of F S(j, k). For every two consecutive claws C i and C i+1 , by Lemma 1, we have either two edges of M joining the external vertices of C i to those of C i+1 or none . We have 3 ways to choose 2 edges between C i and C i+1 , each choice of these two edges can be completed in a unique way in a perfect matching of the subgraph C i ∪ C i+1 . Hence we get easily that the number of perfect matchings of type 2 in F S(j, k)

(j ∈ {1, 2, 3}) is µ 2 (j, k) = 2 × 3 k 2 .
From Theorem 5 and Theorem 6 we deduce :

Corollary 7 The numbers µ(i, k) of perfect matchings of F S(i, k) (i ∈ {1, 2, 3}) are given by: When k is odd

• µ(2, k) = 2 k • µ(1, k) = 2 k + 1 • µ(3, k) = 2 k -2
When k is even

• µ(2, k) = 2 × 3 k 2 + 2 k • µ(1, k) = 2 × 3 k 2 + 2 k -1 • µ(3, k) = 2 × 3 k 2 + 2 k + 2
3 Some structural results about perfect matchings of FS(j, k)

3.1 Perfect matchings of type 1

Lemma 8 Let M be a perfect matching of type 1 of G = F S(j, k). Then the 2-factor G \ M has exactly one or two cycles and each cycle of G \ M has at least one vertex in each claw C i (i ∈ Z k ).

Proof Let M be a perfect matching of type 1 in G. Let us consider the claw

C i for some i in Z k . Assume without loss of generality that the edge of M contained in C i is t i x i . The cycle of G \ M visiting x i comes from C i-1
, crosses C i by using the vertex x i and goes to C i+1 . By Lemma 1, the path y i t i z i is contained in a cycle of G \ M . The two edges incident to y i and z i joining C i to C i-1 (as well as those joining C i to C i+1 ) are not contained both in M (since M has type 1). Thus, the cycle of G \ M containing y i t i z i comes from C i-1 , crosses C i and goes to C i+1 . Thus, we have at most two cycles in G \ M , as claimed, and we can note that each claw must be visited by these cycles.

Denition 9 Let us suppose that M is a perfect matching of type 1 in G = F S(j, k) such that the 2-factor G \ M has exactly two cycles Γ 1 and Γ 2 .

A claw C i intersected by three vertices of Γ 1 (respectively Γ 2 ) is said to be Γ 1 -major (respectively Γ 2 -major).

Lemma 10 Let M be a perfect matching of type 1 of G = F S(j, k) such that the 2-factor G \ M has exactly two cycles. Then, the lengths of these two cycles have the same parity as k, and those lengths are distinct when k is odd.

Proof Let Γ 1 and Γ 2 be the two cycles of G \ M . By Lemma 8, for each i in Z k these two cycles must cross the claw C i . Let k 1 be the number of Γ 1 -major claws and let k 2 be the number of Γ 2 -major claws. We have

k 1 + k 2 = k, l(Γ 1 ) = 3k 1 + k 2 and l(Γ 2 ) = 3k 2 + k 1 .
When k is odd, we must have either k 1 odd and k 2 even, or k 1 even and k 2 odd. Then Γ 1 and Γ 2 have distinct odd lengths. When k is even, we must have either k 1 and k 2 even, or k 2 and k 1 odd. Then Γ 1 and Γ 2 have even lengths.

Lemma 11 Let M be a perfect matching of type 1 of G = F S(j, k) such that the 2-factor G \ M has exactly two cycles Γ 1 and Γ 2 . Suppose that there are two consecutive Γ 1 -major claws C j and C j+1 with j ∈ Z k \ {k -1}. Then there is a perfect matching M ′ of type 1 such that the 2-factor G \ M ′ has exactly two cycles Γ ′ 1 and Γ ′ 2 having the following properties:

a) for i ∈ Z k \ {j, j + 1} C i is Γ ′ 2 -major if and only if C i is Γ 2 -major, b) C j and C j+1 are Γ ′ 2 -major, c) l(Γ ′ 1 ) = l(Γ 1 ) -4 and l(Γ ′ 2 ) = l(Γ 2 ) + 4.
Proof Consider the claws C j and C j+1 . Since C j is a Γ 1 -major claw suppose without loss of generality that t j z j belongs to M and that Γ 1 contains the path

x ′ j-1 x j t j y j y j+1 where x ′ j-1 denotes the neighbour of x j in C j-1 (then x j x j+1 belongs to M ). Since C j+1 is Γ 1 -major and Γ 2 goes through C j and C j+1 , the cycle Γ 1 must contain the path y j+1 t j+1 x j+1 x ′ j+2 where x ′ j+2 denotes the neighbour of x j+1 in C j+2 (then M contains t j+1 z j+1 and y j+1 y ′ j+2 ). Denote by P 1 the path x ′ j-1 x j t j y j y j+1 t j+1 x j+1 x ′ j+2 . Note that Γ 2 contains the path

P 2 = z ′ j-1 z j z j+1 z ′
j+1 where z ′ j-1 and z ′ j+1 are dened similarly. See to the left part of Figure 2.

Let us perform the following local transformation: delete x j x j+1 , t j z j and t j+1 z j+1 from M and add z j z j+1 , t j x j and t j+1 x j+1 . Let M ′ be the resulting perfect matching. Then the subpath P 1 of Γ 1 is replaced by

P ′ 1 = x ′ j-1 x j x j+1 x ′ j+2
and the subpath P 2 of Γ 2 is replaced by P ′ 2 = z ′ j-1 z j t j y j y j+1 t j+1 z j+1 z ′ j+2 (see Figure 2). We obtain a new 2-factor containing two new cycles Γ ′ 1 and Γ ′ 2 . Note that C j and C j+1 are Γ ′ 2 -major claws and for i in

Z k \ {j, j + 1} C i is Γ ′ 2 -major (respectively Γ ′ 1 -major) if and only if C i is Γ 2 -major (respectively Γ 1 -major). The length of Γ 1 (now Γ ′ 1 ) decreases of 4 units while the length of Γ 2 (now Γ ′
2 ) increases of 4 units. The operation depicted in Lemma 11 above will be called a local transformation of type 1.

Lemma 12 Let M be a perfect matching of type 1 of G = F S(j, k) such that the 2-factor G \ M has exactly two cycles Γ 1 and Γ 2 . Suppose that there are three consecutive claws C j , C j+1 and C j+2 with j in Z k \ {k -1, k -2} such that C j and C j+2 are Γ 1 -major and C j+1 is Γ 2 -major. Then there is a perfect matching M ′ of type 1 such that the 2-factor G \ M ′ has exactly two cycles Γ ′ 1 and Γ ′ 2 having the following properties:

a) for i ∈ Z k \ {j, j + 1, j + 2} C i is Γ ′ 2 -major if and only if C i is Γ 2 -major, b) C j and C j+2 are Γ ′ 2 -major and C j+1 is Γ ′ 1 -major, c) l(Γ ′ 1 ) = l(Γ 1 ) -2 and l(Γ ′ 2 ) = l(Γ 2 ) + 2.
Proof Since C j is Γ 1 -major, as in the proof of Lemma 11 suppose that Γ 1 contains the path x ′ j-1 x j t j y j y j+1 (that is edges t j z j and x j x j+1 belong to M ). Since C j+1 is Γ 2 -major the cycle Γ 1 contains the edge y j+1 y j+2 . Then we see that Γ 1 contains the path Q 1 = x ′ j-1 x j t j y j y j+1 y j+2 t j+2 z j+2 z ′ j+3 and that Γ 2 contains the path Q 2 = z ′ j-1 z j z j+1 t j+1 x j+1 x j+2 x ′ j+3 . Note that y j+1 t j+1 , z j+1 z j+2 and t j+2 x j+2 belong to M . Let us perform the following local transformation: delete t j z j , x j x j+1 , z j+1 z j+2 and x j+2 t j+2 from M and add x j t j , z j z j+1 , x j+1 x j+2 and z j+2 t j+2 to M . Let M ′ be the resulting perfect matching. Then the subpath

Q 1 of Γ 1 is replaced by Q ′ 1 = x ′ j-1 x j x j+1 t j+1 z j+1 z j+2 z ′ j+3 and the subpath Q 2 of Γ 2 is replaced by Q ′ 2 = z ′ j-1 z j t j y j y j+1 y j+2 t j+2 x j+2
x ′ j+3 (see Figure 3). We obtain a new 2-factor containing two new cycles named Γ ′ 1 and Γ ′ 2 . Note that C j and C j+2 are now Γ ′ 2major claws and C j+1 is Γ ′ 1 -major. The length of Γ 1 decreases of 2 units while the length of Γ 2 increases of 2 units. It is clear that for The operation depicted in Lemma 12 above will be called a local transformation of type 2.

i ∈ Z k \ {j, j + 1, j + 2} C i is Γ ′ 2 -major (respectively Γ ′ 1 -major) if and only if C i is Γ 2 -major (respec- tively Γ 1 -major).
Lemma 13 Let M be a perfect matching of type 1 of G = F S(j, k) such that the 2-factor G \ M has exactly two cycles Γ 1 and Γ 2 . Suppose that there are three consecutive claws C j , C j+1 and C j+2 with j in Z k \ {k -1, k -2} such that C j+1 and C j+2 are Γ 2 -major and C j is Γ 1 -major. Then there is a perfect matching M ′ of type 1 such that the 2-factor G \ M ′ has exactly two cycles Γ ′ 1 and Γ ′ 2 having the following properties:

a) for i ∈ Z k \ {j, j + 1, j + 2} C i is Γ ′ 2 -major if and only if C i is Γ 2 -major, b) C j and C j+1 are Γ ′ 2 -major and C j+2 is Γ ′
Proof Since C j is Γ 1 -major, as in the proof of Lemma 11 suppose that Γ 1 contains the path x ′ j-1 x j t j y j y j+1 (that is edges t j z j and x j x j+1 belong to M ). Since C j+1 and C j+2 are Γ 2 -major, the unique vertex of C j+1 (respectively C j+2 ) contained in Γ 1 is y j+1 (respectively y j+2 ). Note that the perfect matching M contains the edges t j z j , x j x j+1 , t j+1 y j+1 , z j+1 z j+2 and t j+2 y j+2 . Then the path R 1 = x ′ j-1 x j t j y j y j+1 y j+2 y ′ j+3 is a subpath of Γ 1 and the path R 2 = z ′ j-1 z j z j+1 t j+1 x j+1 x j+2 t j+2 z j+2 z ′ j+3 is a subpath of Γ 2 . See to the left part of Figure 4.

Let us perform the following local transformation: delete t j z j , x j x j+1 , t j+1 y j+1 , z j+1 z j+2 and t j+2 y j+2 from M and add x j t j , z j z j+1 , t j+1 x j+1 , y j+1 y j+2 and t j+2 z j+2 . Let M ′ be the resulting perfect matching. Then the subpath R

1 of Γ 1 is replaced by R ′ 1 = x ′ j-1 x j x j+1 x j+2 t j+2 y j+2 y ′ j+3 and the subpath R 2 of Γ 2 is replaced by R ′ 2 = z ′ j-1 z j t j y j y j+1 t j+1 z j+1 z j+2 z ′ j+3 .
We obtain a new 2factor containing two new cycles named Γ ′ 1 and Γ ′ 2 such that l(Γ ′ 1 ) = l(Γ 1 ) and l(Γ ′ 2 ) = l(Γ 2 ) (see Figure 4). It is clear that for i ∈ Z k \ {j, j + 1, j + 2} C i is Γ ′ 2 -major (respectively Γ ′ 1 -major) if and only if C i is Γ 2 -major (respectively Γ 1 -major). Note that C j and C j+1 are Γ ′ 2 -major and C j+2 is Γ ′ 1 -major. The operation depicted in Lemma 13 above will be called a local transformation of type 3. Lemma 14 Let M be a perfect matching of type 1 of G = F S(j, k) such that the 2-factor G \ M has exactly two cycles Γ 1 and Γ 2 such that l(Γ 1 ) ≤ l(Γ 2 ) and l(Γ 2 ) is as great as possible. Then there exists at most one Γ 1 -major claw.

Proof Suppose, for the sake of contradiction, that there exist at least two Γ 1 -major claws. Since l(Γ 2 ) is maximum, by Lemma 11 these claws are not consecutive. Then consider two Γ 1 -major claws C i and C i+h+1 (with h ≥ 1) such that the h consecutive claws (C i+1 , . . . , C i+h ) are Γ 2 -major. Since l(Γ 2 ) is maximum, by Lemma 12 the number h is at least 2. Then by applying r = ⌊ h 2 ⌋ consecutive local transformations of type 3 (Lemma 13) we obtain a perfect matching M (r) such that the 2-factor G \ M (r) has exactly two cycles Γ 2 ) is maximum, we can conclude by Lemma 11 and by Lemma 12 that h is neither even nor odd, a contradiction.

Perfect matchings of type 2

We give here a structural result about perfect matchings of type 2 in G = F S(j, k).

Lemma 15 Let M be a perfect matching of type 2 of G = F S(j, k) (with k ≥ 4). Then the 2-factor G \ M has exactly one cycle of even length l ≥ k and a set of p cycles of length 6 where l + 6p = 4k (with 0 ≤ p ≤ k 2 ).

Proof Let M be a perfect matching of type 2 in G. By Lemma 1 the number k of claws is even. Let i in Z k such that there are two edges of M between C i-1 and C i . There are no edges of M between C i and C i+1 and two edges of M between C i+1 and C i+2 . We may consider that 0 ≤ i < k -1. For j ∈ {i, i + 2, i + 4, . . .} we denote by e j the unique edge of G \ M having one end vertex in C j-1 and the other in C j . Let us denote by A the set {e i , e i+2 , e i+4 , . . .}. We note that

| A |= k 2 .
Assume without loss of generality that the two edges of M between C i-1 and C i have end vertices in C i which are x i and y i (then z i is the end vertex of e i in C i ). Two cases may now occur.

Case 1: The end vertices in C i+1 of the two edges of M between C i+1 and C i+2 are x i+1 and y i+1 (then z i+1 is the end vertex of e i+2 in C i+1 ). In that case the 2factor G\M contains the cycle of length 6 x i x i+1 t i+1 y i+1 y i t i while the edge z i z i+1 of G \ M connects e i and e i+2 .

Case 2: The end vertices in C i+1 of the two edges of M between C i+1 and C i+2 are y i+1 and z i+1 (respectively x i+1 and z i+1 ). Then x i+1 (respectively y i+1 ) is the end vertex of e i+2 in C i+1 . In that case the edges e i and e i+2 are connected in G \ M by the path z i z i+1 t i+1 y i+1 y i t i x i x i+1 (respectively z i z i+1 t i+1 x i+1 x i t i y i y i+1 ).

The same reasoning can be done for {e i+2 , e i+4 }, {e i+4 , e i+6 }, and so on. Then, we see that the set A is contained in a unique cycle Γ of G \ M which crosses each claw. Thus, the length l of Γ is at least k. More precisely, each e j in A contributes for 1 in l, in Case 1 the edge z i z i+1 contributes for 1 in l and in Case 2 the path z i z i+1 t i+1 y i+1 y i t i x i x i+1 contributes for 7 in l. Let us suppose that Case 1 appears p times (0 ≤ p ≤ k 2 ), that is to say G \ M contains p cycles of length 6. Since Case 2 appears k 2 -p times, the length of Γ is

l = k 2 + p + 7( k 2 -p) = 4k -6p.
Remark 16 If k is even then by Lemmas 8, 10 and 15 F S(j, k) has an even 2-factor. That is to say F S(j, k) is a cubic 3-edge colourable graph .

4 Perfect matchings and hamiltonian cycles of F (j, k)

Perfect matchings of type 1 and hamiltonicity

Theorem 17 Let M be a perfect matching of type 1 of G = F S(j, k). Then the 2-factor G \ M is a hamiltonian cycle except for k odd and j = 2, and for k even and j = 1 or 3.

Proof Suppose that there exists a perfect matching M of type 1 of G such that G \ M is not a hamiltonian cycle. By Lemma 8 and Lemma 10 the 2factor G \ M is made of exactly two cycles Γ 1 and Γ 2 whose lengths have the same parity as k. Without loss of generality we suppose that l(Γ 1 ) ≤ l(Γ 2 ).

Assume moreover that among the perfect matchings of type 1 of G such that the 2-factor G \ M is composed of two cycles, M has been chosen in such a way that the length of the longest cycle Γ 2 is as great as possible. By Lemma 14 there exists at most one Γ 1 -major claw.

Case 1: There exists one Γ 1 -major claw.

Without loss of generality, suppose that C 0 is intersected by Γ 1 into {y 0 , t 0 , x 0 } and that y ′ k-1 y 0 belongs to Γ 1 . Since for every i = 0 the claw C i is Γ 2 -major, Γ 1 contains the vertices y 0 , t 0 , x 0 , x 1 , x 2 , . . . , x k-1 which means that x ′ k-1 = y 0 .

• If k = 2r + 1 with r ≥ 1 then Γ 2 contains the path z 0 z 1 t 1 y 1 y 2 t 2 z 2 ...z 2r-1 t 2r-1 y 2r-1 y 2r t 2r z 2r . Thus, y 0 x k-1 , x 0 y k-1 , z 0 z k-1 are edges of G. This means that ∪ i=k-1 i=0 {C i \ {t i }} induces two cycles, that is to say j = 2 and G = F S(2, k). • If k = 2r + 2 with r ≥ 1 then Γ 2 contains the path z 0 z 1 t 1 y 1 y 2 t 2 z 2 ...z 2r-1 t 2r-1 y 2r-1 y 2r t 2r z 2r z 2r+1 t 2r+1 y 2r+1 .
Thus, x 0 z k-1 , y 0 x k-1 and z 0 y k-1 are edges. This means that ∪ i=k-1 i=0 {C i \ {t i }} induces one cycle, that is to say j = 1 and G = F S(1, k).

When j = 1 or 3 and k is odd by Theorem 17 F S(j, k) is hamiltonian. If k is even then by Lemmas 8, 10 and 15 F S(j, k) has an even 2-factor.

4.2

Perfect matchings of type 2 and hamiltonicity At this point of the discourse one may ask what happens for perfect matchings of type 2 in F S(j, k) (k even). Can we characterize and count perfect matchings of type 2, complementary 2-factor of which is a hamiltonian cycle ? An armative answer shall be given.

Let us consider a perfect matching M of type 2 in F S(j, 2p) with p ≥ 2.

Suppose that there are no edges of M between C 2i-1 and C 2i (for any i ≥ 1), that is to say M is a matching of type 2.0 (see Denition 2). Consider two consecutive claws C 2i and C 2i+1 (0 ≤ i ≤ p -1). There are three cases:

Case (x):

{y 2i y 2i+1 , z 2i z 2i+1 } ⊂ M (then, M ∩(C 2i ∪C 2i+1 ) = {x 2i t 2i , x 2i+1 t 2i+1 }).
Case (y):

{x 2i x 2i+1 , z 2i z 2i+1 } ⊂ M (then, M ∩(C 2i ∪C 2i+1 ) = {y 2i t 2i , y 2i+1 t 2i+1 }).
Case (z):

{x 2i x 2i+1 , y 2i y 2i+1 } ⊂ M (then, M ∩(C 2i ∪C 2i+1 ) = {z 2i t 2i , z 2i+1 t 2i+1 }).
The subgraph induced on C 2i ∪C 2i+1 is called a block. In Case (x) (respectively Case (y), Case (z)) a block is called a block of type X (respectively block of type Y , block of type Z). Then F S(j, 2p) with a perfect matching M of type 2.0 can be seen as a sequence of p blocks properly connected. In other words, a perfect matching M of type 2 in F S(j, 2p) is entirely described by a word of length p on the alphabet of three letters {X, Y, Z}. The block C 0 ∪ C 1 is called initial block and the block C 2p-1 ∪ C 2p is called terminal block. These extremal blocks are not considered here as consecutive blocks.

By Lemma 15, F S(j, 2p) \ M has no 6-cycles if and only if F S(j, 2p) \ M is a unique even cycle. It is an easy matter to prove that two consecutive blocks do not induce a 6-cycle if and only if they are not of the same type. Then the possible congurations for two consecutive blocks are XY , XZ, Y X, Y Z, ZX and ZY . To eliminate a possible 6-cycle in C 0 ∪ C 2p-1 we have to determine for every j ∈ {1, 2, 3} the forbidden extremal congurations. An extremal conguration shall be denoted by a word on two letters in {X, Y, Z} such that the left letter denotes the type of the initial block C 0 ∪ C 1 and the right letter denotes the type of the terminal block C 2p-1 ∪ C 2p . We suppose that the extremal blocks are connected for j = 1 by the edges x 2p-1 z 0 , y 2p-1 x 0 and z 2p-1 y 0 , for j = 2 by the edges x 2p-1 x 0 , y 2p-1 z 0 and z 2p-1 y 0 and for j = 3 by the edges x 2p-1 x 0 , y 2p-1 y 0 and z 2p-1 z 0 . Then, it is easy to verify that we have the following result.

Lemma 21 Let M be a perfect matching of type Theorem 22 The numbers µ ′ 2 (j, 2p) of perfect matchings of type 2 complementary to hamiltonian cycles in F S(j, 2p) (j ∈ {1, 2, 3}) are given by: Proof Consider, as previously, perfect matchings of type if α = β, and a 2 αβ = 1 if α = β. If α = β we have to solve the recurrent sequence : u 2 = 2 and u p = 2 p-1 -u p-1 for p ≥ 3. If α = β we have to solve y 0 y 1 belong to M B . Since M B is a strong matching, t 0 z 0 and t 1 z 1 belong to M \ M B = M R . This is impossible because M R is also a strong matching. By symmetry there are no two edges of M R between C 0 and C 1 . Then there is one edge of M B between C 0 and C 1 , x 0 x 1 for instance, and one edge of M R between C 0 and C 1 , y 0 y 1 for instance. Since M B and M R are strong matchings, there is no edge of M in C 0 ∪ C 1 , a contradiction. Thus, M is a perfect matching of type 1.

µ ′ 2 (1, 2p) = 2 p+1 + (-1) p+1 2,
Lemma 25 If G = F S(j, k) is a Jaeger's graph (with k ≥ 3) then either (j = 1 and k ≡ 1 or 2 (mod 3)) or (j = 3 and k ≡ 0 (mod 3)).

Proof Let M = M B ∪ M R be a Jaeger's matching of G. By Lemma 24 M is a perfect matching of type 1. Suppose without loss of generality that M B ∩ E(C 0 ) = {x 0 t 0 }. Since M B is a strong matching there is no edge of M B between C 0 and C 1 . Suppose, without loss of generality, that the edge in M R joining C 0 to C 1 is y 0 y 1 . Consider the claws C 0 , C 1 and C 2 . Since M B and M R are strong matchings, we can see that the choices of x 0 t 0 ∈ M B and y 0 y 1 ∈ M R xes the positions of the other edges of M B and M R . More precisely, {t 1 z 1 , y 2 t 2 } ⊂ M B and {x 1 x 2 , z 2 z ′ 3 } ⊂ M R . This unique conguration is depicted in Figure 5. If k ≥ 4 then we see that z 2 z 3 ∈ M R , x 3 t 3 ∈ M B , and y 3 y ′ 4 ∈ M R . So, the local situation in C 3 is similar to that in C 0 , and we can see that there is a unique Jaeger's matching M = M B ∪ M R such that x 0 t 0 ∈ M B and y 0 y 1 ∈ M R in the graph F S(j, k). We have to verify the coherence of the connections between the claws C k-1 and C 0 . We note that

M B = M ∩ (∪ i=k-1 i=0 E(C i )) and M R is a strong matching included in the 2-factor induced by ∪ i=k-1 i=0 {V (C i ) \ {t i }}. Case 1: k = 3p with p ≥ 1. We have x 0 t 0 ∈ M B , y k-1 t k-1 ∈ M B , x k-2 x k-1 ∈ M R and z ′ k-1 z 0 = z k-1 z ′ 0 ∈ M R (that is, z k-1 z 0 ∈ M R )
. Thus, z k-1 z 0 , y k-1 y 0 and x k-1 x 0 are edges of F S(j, 3p) and we must have j = 3.

Case 2: k = 3p + 1 with p ≥ 1.

and such a pair {e, f } corresponds exactly to one Jaeger's matching. Then, a Jaeger's graph F S(j, k) (with k ≥ 3) has exactly 6 Jaeger's matchings.

Remark 28 The Berge-Fulkerson Conjecture states that if G is a bridgeless cubic graph, then there exist six perfect matchings M 1 , . . . , M 6 of G (not necessarily distinct) with the property that every edge of G is contained in exactly two of M 1 , . . . , M 6 (this conjecture is attributed to Berge in [START_REF] Seymour | On multi-colourings of cubic graphs, and conjectures of Fulkerson and Tutte[END_REF] but appears in [START_REF] Fulkerson | Blocking and anti-blocking pairs of polyhedra[END_REF]). Using each colour of a cubic 3-edge colourable graph twice, we see that such a graph fulls the Berge-Fulkerson Conjecture. Very few is known about this conjecture excepted that it holds for the Petersen graph and for cubic 3-edge colourable graphs. So, Berge-Fulkerson Conjecture holds for Jaeger's graphs, but generally we do not know if we can nd six distinct perfect matchings. We remark that if F S(j, k), with k ≥ 3, is a Jaeger's graph then its six Jaeger's matchings are such that every edge is contained in exactly two of them.

2-factor hamiltonian cubic graphs

Recall that a simple graph of maximum degree d > 1 with edge chromatic number equal to d is said to be a Class 1 graph. For any d-regular simple graph (with d > 1) of even order and of Class 1, for any minimum edgecolouring of such a graph, the set of edges having a given colour is a perfect matching (or 1-factor). Such a regular graph is also called a 1-factorable graph.

A Class 1 d-regular graph of even order is strongly hamiltonian or perfectly 1factorable (or is a Hamilton graph in the Kotzig's terminology [START_REF] Kotzig | Balanced colourings and the four colour conjecture[END_REF]) if it has an edge colouring such that the union of any two colours is a hamiltonian cycle. Such an edge colouring is said to be a Hamilton decomposition in the Kotzig's terminology. In [START_REF] Kotzig | Construction for hamiltonian graphs of degree three (in russian)[END_REF] by using two operations ρ and π (described also in [START_REF] Kotzig | Balanced colourings and the four colour conjecture[END_REF]) and starting from the θ-graph (two vertices joined by three parallel edges) he obtains all strongly hamiltonian cubic graphs, but these operations do not always preserve planarity. In his paper [START_REF] Kotzig | Balanced colourings and the four colour conjecture[END_REF] he describes a method for constructing planar strongly hamiltonian cubic graphs and he deals with the relation between strongly hamiltonian cubic graphs and 4-regular graphs which can be decomposed into two hamiltonian cycles. See also [START_REF] Kotzig | Quelques problèmes ouverts concernant les graphes fortement hamiltoniens[END_REF] and a recent work on strongly hamiltonian cubic graphs [START_REF] Bonvicini | On perfectly one-factorable cubic graphs[END_REF] in which the authors give a new construction of strongly hamiltonian graphs.

A Class 1 regular graph such that every edge colouring is a Hamilton decomposition is called a pure Hamilton graph by Kotzig [START_REF] Kotzig | Balanced colourings and the four colour conjecture[END_REF]. Note that K 4 is a pure Hamilton graph and every cubic graph obtained from K 4 by a sequence of triangular extensions is also a pure Hamilton cubic graph. In the paper [START_REF] Kotzig | Balanced colourings and the four colour conjecture[END_REF] of Kotzig, a consequence of his Theorem 9 (p.77) concerning pure Hamilton graphs is that the family of pure Hamilton graphs that he exhibits is precisely the family obtained from K 4 by triangular extensions. Are there others pure Hamilton cubic graphs ? The answer is "yes".

We remark that 2-factor hamiltonian cubic graphs dened above (see Denition 18) are pure Hamilton graphs (in the Kotzig's sense) but the converse is false because K 4 is 2-factor hamiltonian and the pure Hamilton cubic graph on 6 vertices obtained from K 4 by a triangular extension (denoted by P R 3 ) is not 2-factor hamiltonian. Observe that the operation of triangular extension preserves the property "pure Hamilton", but does not preserve the property "2-factor hamiltonian". The Heawood graph H 0 (on 14 vertices) is pure Hamiltonian, more precisely it is 2-factor hamiltonian (see [START_REF] Funk | On minimally one-factorable r-regular bipartite graphs[END_REF] Proposition 1.1 and Remark 2.7). Then, the graphs obtained from the Heawood graph H 0 by triangular extensions are also pure Hamilton graphs.

A minimally 1-factorable graph G is dened by Labbate and Funk [START_REF] Funk | On minimally one-factorable r-regular bipartite graphs[END_REF] as a Class 1 regular graph of even order such that every perfect matching of G is contained in exactly one 1-factorization of G. In their article they study bipartite minimally 1-factorable graphs and prove that such a graph G has necessarily a degree d ≤ 3. If G is a minimally 1-factorable cubic graph then the complementary 2-factor of any perfect matching has a unique decomposition into two perfect matchings, therefore this 2-factor is a hamiltonian cycle of G, that is G is 2-factor hamiltonian. Conversely it is easy to see that any 2-factor hamiltonian cubic graph is minimally 1-factorable. The complete bipartite graph K 3,3 and the Heawood graph H 0 are examples of 2-factor hamiltonian bipartite graph given by Labbate and Funk. Starting from H 0 , from K 1,3 and from three copies of any tree of maximum degree 3 and using three operations called amalgamations the authors exhibit an innite family of bipartite 2-factor hamiltonian cubic graphs, namely the poly -HB -R -R 2 graphs (see [START_REF] Funk | On minimally one-factorable r-regular bipartite graphs[END_REF] for more details). Except H 0 , these graphs are exactly cyclically 3-edge connected. Others structural results about 2-factor hamiltonian bipartite cubic graph are obtained in [START_REF] Labbate | On 3-cut reductions of minimally 1-factorable cubic bigraphs[END_REF], [START_REF] Labbate | Characterizing minimally 1-factorable r-regular bipartite graphs[END_REF]. These results have been completed and a simple method to generate 2-factor hamiltonian bipartite cubic graphs was given in [START_REF] Funk | 2-factor hamiltonian graphs[END_REF].

Proposition 29 (Lemma 3.3, [START_REF] Funk | 2-factor hamiltonian graphs[END_REF]) Let G be a 2-factor hamiltonian bipartite cubic graph. Then G is 3-connected and | V (G) |≡ 2 (mod 4).

Let G 1 and G 2 be disjoint cubic graphs, x ∈ v(G 1 ), y ∈ v(G 2 ). Let x 1 , x 2 , x 3 (respectively y 1 , y 2 , y 3 ) be the neighbours of x in G 1 (respectively, of y in G 2 ). The cubic graph G such that V (G) = (V (G 1 )\{x})∪(V (G 2 )\{y}) and E(G) = (E(G 1 )\{x 1 x, x 2 x, x 3 x})∪(E(G 2 )\{y 1 y, y 2 y, y 3 y})∪{x 1 y 1 , x 2 y 2 , x 3 y 3 } is said to be a star product and G is denoted by (G 1 , x) * (G 2 , y). Since {x 1 y 1 , x 2 y 2 , x 3 y 3 } is a cyclic edge-cut of G, a star product of two 3-connected cubic graphs has cyclic edge-connectivity 3.

Fig. 1 .

 1 Fig. 1. Four consecutive claws

Fig. 2 .

 2 Fig. 2. Local transformation of type 1

Fig. 3 .

 3 Fig. 3. Local transformation of type 2

Fig. 4 .

 4 Fig. 4. Local transformation of type 3

  l(Γ 1 ) and l(Γ (r)2 ) = l(Γ 2 ) and such that C i+2⌊ h 2 ⌋ and C i+h+1 are Γ

µ ′ 2 ( 2 ,

 22 2p) = 2 p+1 , and µ ′ 2 (3, 2p) = 2 p+1 + (-1) p 4.

Fig. 5 .

 5 Fig. 5. Strong matchings M B (bold edges) and M R (dashed edges)

5

 5

  2.0 of G = F S(j, 2p) (with p ≥ 2) such that the 2-factor G\M is a hamiltonian cycle. Then the forbidden extremal congurations are XY , Y Z and ZX for F S(1, 2p), XX, Y Z and ZY for F S(2, 2p), and XX, Y Y and ZZ for F S(3, 2p).Thus, any perfect matching M of type 2.0 of F S(j, 2p) such that the 2-factor G \ M is a hamiltonian cycle is totally characterized by a word of length p on the alphabet {X, Y, Z} having no two identical consecutive letters and such that the sub-word [initial letter][terminal letter] is not a forbidden conguration. Then, we are in position to obtain the number of such perfect matchings in F S(j, 2p). Let us denote by µ ′ 2.0 (j, 2p) (respectively µ ′ the number of perfect matchings of type 2.0 (respectively type 2.1, type 2) complementary to a hamiltonian cycle in F S(j, 2p). Clearly µ ′

	2.1 (j, 2p), 2 (j, 2p)) 2 (j, 2p) = µ ′ µ ′ 2.0 (j, 2p) + µ ′ 2.1 (j, 2p) and µ ′ 2.0 (j, 2p) = µ ′ 2.1 (j, 2p).

  2.0. Let α and β be two letters in {X, Y, Z} (not necessarily distinct). Let A p αβ be the set of words of length p on {X, Y, Z} having no two consecutive identical letters, beginning by α and ending by a letter distinct from β. Denote the number of words in A p αβ by a p αβ . Let B p αβ be the set of words of length p on {X, Y, Z} having no two consecutive identical letters, beginning by α and ending by β. Denote by

	b p αβ the number of words in B p αβ .
	Clearly, the number of words of length p having no two consecutive identical letters and beginning by α is 2 p-1 . Then a p αβ + b p αβ = 2 p-1 . The deletion of the
	last β of a word in B p αβ gives a word in A p-1 αβ and the addition of β to the right
	of a word in A p-1 αβ gives a word in B p αβ .
	Thus b p αβ = a p-1 αβ and for every p ≥ 3 a p αβ = 2 p-1 -a p-1 αβ . We note that a 2 αβ = 2

Preprint submitted to Discussiones Math. -Graph Theory

-major, c) l(Γ ′ 1 ) = l(Γ 1 ) and l(Γ ′

) = l(Γ 2 ).

The authors specify that a smallest counterexample to Conjecture 31 is a cyclically 4-edge connected cubic graph of girth at least 6, and that to show this result it would suce to prove that H 0 is the only 2-factor hamiltonian cyclically 4-edge connected bipartite cubic graph of girth at least 6. Note that some results have been generalized in [1].

To conclude, we may ask what happens for non bipartite 2-factor hamiltonian cubic graphs. Recall that K 4 and F S(1, 3) (the "Triplex Graph" of Robertson, Seymour and Thomas [15]) are 2-factor hamiltonian cubic graphs. By Corollary 19 the graphs F S(j, k) with k odd and j = 1 or 3 introduced in this paper form a new innite family of non bipartite 2-factor hamiltonian cubic graphs. We remark that they are cyclically 6-edge connected. Can we generate others families of non bipartite 2-factor hamiltonian cubic graphs ? Since P R 3 (the cubic graph on 6 vertices obtained from K 4 by a triangular extension) is not 2-factor hamiltonian and P R 3 = K 4 * K 4 , the star product operation is surely not a possible tool.

Case 2: There is no Γ 1 -major claw.

Suppose that x 0 belongs to Γ 1 . Then, Γ 1 contains x 0 , x 1 , ..., x k-1 .

• If k = 2r + 1 with r ≥ 1 then Γ 2 contains the path y 0 t 0 z 0 z 1 t 1 y 1 y 2 ...z 2r-1 t 2r-1 y 2r-1 y 2r t 2r z 2r . Thus, x 0 x k-1 , y 0 z k-1 and z 0 y k-1 are edges of G and the set ∪ i=k-1 i=0 {C i \ {t i }} induces two cycles,that is to say j = 2 and G = F S(2, k).

• If k = 2r + 2 with r ≥ 1 then Γ 2 contains the path y 0 t 0 z 0 z 1 t 1 y 1 y 2 ...y 2r t 2r z 2r z 2r+1 t 2r+1 y 2r+1 . Thus, x 0 x k-1 , y 0 y k-1 and z 0 z k-1 are edges. This means that ∪ i=k-1 i=0 {C i \ {t i }} induces three cycles, that is to say j = 3 and G = F S(3, k).

Denition 18 A cubic graph G is said to be 2-factor hamiltonian [START_REF] Funk | 2-factor hamiltonian graphs[END_REF] if every 2-factor of G is a hamiltonian cycle (or equivalently, if for every perfect matching M of G the 2-factor G \ M is a hamiltonian cycle).

By Theorem 17 for any odd k ≥ 3 and j ∈ {1, 3} or for any even k and j = 2, and for every perfect matching M of type 1 in F S(j, k) the 2-factor F S(j, k) \ M is a hamiltonian cycle. By Lemma 15 F S(2, k) (k ≥ 4) may have a perfect matching M of type 2 such that the 2-factor F S(2, k) \ M is not a hamiltonian cycle (it may contains cycles of length 6).

Then we have the following.

Corollary 19 A graph G = F S(j, k) is 2-factor hamiltonian if and only if k is odd and j = 1 or 3.

We note that F S(1, 3) is the "Triplex Graph" of Robertson, Seymour and Thomas [START_REF] Robertson | Excluded minor in cubic graphs[END_REF]. We shall examine other known results about 2-factor hamiltonian cubic graphs in Section 5.

Corollary 20 The chromatic index of a graph G = F S(j, k) is 4 if and only if j = 2 and k is odd.

Proof When j = 2 and k is odd, any 2-factor must have at least two cycles, by Theorem 17. Then Lemma 10 implies that any 2-factor is composed of two odd cycles. Hence G has chromatic index 4.

the recurrent sequence : v 2 = 1 and v p = 2 p-1 -v p-1 for p ≥ 3. Then we obtain u p = 2 3 (2 p-1 + (-1) p ) and v p = 1 3 (2 p + (-1) p+1 ) for p ≥ 2.

By Lemma 21

) we obtain the announced results.

Remark 23 We see that µ ′ 2 (j, 2p) ≃ 2 p+1 and this is to compare with the number µ 2 (j, 2p) = 2 × 3 p of perfect matchings of type 2 in F S(j, 2p) (see backward in Section 2).

Strong matchings and Jaeger's graphs

For a given graph G = (V, E) a strong matching (or induced matching) is a matching S such that no two edges of S are joined by an edge of G. That is, S is the set of edges of the subgraph of G induced by the set V (S). We consider cubic graphs having a perfect matching which is the union of two strong matchings that we call Jaeger's graph (in his thesis [START_REF] Jaeger | Etude de quelques invariants et problèmes d'existence en théorie de graphes[END_REF] Jaeger called these cubic graphs equitable). We call Jaeger's matching a perfect matching M of a cubic graph G which is the union of two strong matchings M B and M R . Set B = V (M B ) (the blue vertices) and R = V (M R ) (the red vertices).

An edge of G is said mixed if its end vertices have distinct colours. Since the set of mixed edges is E(G) \ M , the 2-factor G \ M is even and

Thus, every Jaeger's graph G is a cubic 3-edge colourable graph and for any Jaeger's matching

See, for instance, [3] and [START_REF] Fouquet | On isomorphic linear partition in cubic graphs[END_REF] for some properties of these graphs.

In this subsection we determine the values of j and k for which a graph F S(j, k)

is a Jaeger's graph.

Lemma 24 If G = F S(j, k) is a Jaeger's graph (with k ≥ 3) and M = M B ∪ M R is a Jaeger's matching of G then M is a perfect matching of type 1.

Proof Suppose that M is of type 2 and suppose without loss of generality that there are two edges of M between C 0 and C 1 , for instance x 0 x 1 and y 0 y

We have

. Thus, y k-1 z 0 , x k-1 y 0 and z k-1 x 0 are edges of F S(j, 3p + 1) and we must have j = 1.

Case 3: k = 3p + 2 with p ≥ 1.

We have

x 0 and z k-1 y 0 are edges of F S(j, 3p + 2) and we must have j = 1.

Remark 26 It follows from Lemma 25 that for every k ≥ 3 the graph F S(2, k) is not a Jaeger's graph. This is obvious when k is odd, since the ower snarks have chromatic index 4.

Then, we obtain the following.

Theorem 27 For j ∈ {1, 2, 3} and k ≥ 2, the graph G = F S(j, k) is a Jaeger's graph if and only if either k ≡ 1 or 2 (mod 3) and j = 1, or k ≡ 0 (mod 3) and j = 3.

Moreover, F S(1, 2) has 3 Jaeger's matchings and for k ≥ 3 a Jaeger's graph G = F S(j, k) has exactly 6 Jaeger's matchings.

Proof For k = 2 we remark that F S(1, 2) (that is the cube) has exactly three distinct Jaeger's matchings M 1 , M 2 and M 3 . Following our notations: M 1 = {x 0 t 0 , t 1 z 1 } ∪ {y 0 y 1 , z 0 x 1 }, M 2 = {z 0 t 0 , t 1 y 1 } ∪ {y 0 z 1 , x 0 x 1 } and M 3 = {y 0 t 0 , t 1 x 1 } ∪ {z 0 z 1 , x 0 y 1 }.

For k ≥ 3, by Lemma 25, condition ( * ) (j = 1 and k ≡ 1 or 2 (mod 3)) or (j = 3 and k ≡ 0 (mod 3))

is a necessary condition for F S(j, k) to be a Jaeger's graph.

Consider the function

Dene similarly Φ X,Z and Φ Y,Z . For j = 1 or 3 these functions are automorphisms of F S(j, k). Thus, the process described in the proof of Lemma 25 is a constructive process of all Jaeger's matchings in a graph F S(j, k) (with k ≥ 3) verifying condition (*).

We remark that for any choice of an edge e of C 0 to be in M B there are two distinct possible choices for an edge f between C 0 and C 1 to be in M R , Proposition 30 (Proposition 3. Then, taking iterated star products of K 3,3 and the Heawood graph H 0 an innite family of 2-factor hamiltonian cubic graphs is obtained. These graphs (excepted K 3,3 and H 0 ) are exactly cyclically 3-edge connected. In [START_REF] Funk | 2-factor hamiltonian graphs[END_REF] the authors conjecture that the process is complete.

Conjecture ! (Funk, Jackson, Labbate, Sheehan (2003) [START_REF] Funk | 2-factor hamiltonian graphs[END_REF]) Let G be a bipartite 2-factor hamiltonian cubic graph. Then G can be obtained from K 3,3 and the Heawood graph H 0 by repeated star products.


Tan Viet

Anh Truong

Gautier Hattenberger
email: gautier.hattenberger@enac.fr

Catherine Ronfle-Nadaud
email: catherine.ronfle-nadaud@enac.fr

The cooperation between Unmanned Aerial Vehicles using a mission planner

The goal of the paper is to perform a cooperation between multiple Unmanned Aerial Vehicles via a mission planner. This mission planner is a generic one which could be used not only for UAV domain but also for other human-machine systems such as maritime, automotive and manned aerial vehicle (MAV). It is a completed mission planner including the high level task planning and low level of collision avoidance. Our first scenario shows the cooperation between UAVs during doing tasks when the number of UAVs is greater than that of tasks. The later scenarios will be the cooperation between two types of UAV: fixed-wing aircraft and rotor-craft for different types of task. This mission planner contributes a new module to both simulation platform and real demonstrator of relevant projects.

I. INTRODUCTION

In UAV domain, the mission planning issue has been studied longtime ago. They focused not only on trajectory planning [START_REF] Gancet | Task planning and control for a multi-UAV system: architecture and algorithms[END_REF], [START_REF] Mosteo | Multi-robot routing under limited communication range[END_REF] but also on task planning [START_REF] Bellingham | Multi-task allocation and path planning for cooperating UAVs[END_REF], [START_REF] Lagoudakis | Simple auctions with performance guarantees for multi-robot task allocation[END_REF]. For one agent, we plan to do all the sequence of task at the desired time [START_REF] Chanthery | Planning algorithms for autonomous aerial vehicle[END_REF]. The planning is more complicated when we have many tasks to do, then multi agents have to be considered [START_REF] Gerkey | A formal analysis and taxonomy of task allocation in multi-robot systems[END_REF]. Recently, the online task planning [START_REF] Pongpunwattana | Real-time planning for multiple autonomous vehicles in dynamic uncertain environments[END_REF] is also conducted after dealing with the offline one [START_REF] Grøtli | Task assignment for cooperating UAVs under radio propagation path loss constraints[END_REF]. The goal of this paper is to present an online mission planner for multi-UAVs. This mission planner contributes a module to Paparazzi platform.

Paparazzi is an open source platform used to control and simulate offline/online one/multi UAV for different types of UAVs as well as tasks [START_REF] Brisset | The paparazzi solution[END_REF], [11]. So with this new module, the Paparazzi users could perform different kinds of the cooperation between multiple UAVs.

The cooperation of multiple UAVs has been studied few years ago. Among different problems, the mission planning is merged as the main activity. For this problem, we resolve how to allocate available tasks to available aircrafts at one moment of time. Many algorithms have been proposed for mathematical and applied study. In this paper, we propose our mission planner developed by ENAC. This mission planner is used for planning the cooperation between multiple UAVs. The UAV can be the fixed-wing aircraft or rotor-craft depending on the goal type. All the simulation were done by using the Paparazzi simulation platform. The mission planner is integrated as a module in this platform. In the first simulation, we dealt with the scenario of three aircrafts doing two surveillance tasks.

All our scenarios aim the cooperation between UAVs, the autonomy of themselves, and the efficiency of distributing UAVs to tasks.

The paper is organized as following. We first describe our mission planner in section II. The demonstrator architecture and scenario are presented in section III and IV respectively. Finally, we show the simulation results in section V.

II. MISSION PLANNER

A. Mission formulation

We define the mission as the following data structures:

Mission = {Waypoint,

B. High level task planning

Let's consider that we have a set of aircrafts and a set of tasks. The insertion of tasks to aircrafts is based on the time to go from the current position of the aircrafts to the position of tasks and the priority of these tasks. So the multiple criteria function for task allocation is described as the following:

J j(j=1..N ) = min i=1..M ( (t ij , p j )) (1)
where N and M are the number of tasks and aircrafts respectively, p j is the priority factor of the task j and the t ij is the travel time to go from the position of the aircraft i to the position of the task j.

The mission planner algorithm is shown by Algorithm 3 in which we use task planning Algorithm 1 and task replanning Algorithm 2. Our current mission planner works in the case where the number of tasks N is smaller than the number of aircrafts M . The generic one will be studied in the scope of this research.

C. Low level planning of collision avoidance

The "TCAS" (Traffic Collision Avoidance System) technique [START_REF] Brisset | Multi-UAV Control with the Paparazzi System[END_REF], [11] is used for low level of mission planner including the collision avoidance between UAVs. As the first Step 1 Parser XML misision file to mission data structure;

Step 2 Initialization: Sending all aircrafts to their standby points; step, we used the vertical collision avoidance module. It means that when there is a conflict between two UAVs, they will change their altitudes to guarantee their desired trajectories.

Another collision avoidance method such as horizontal collision avoidance (based on the potential force field method) will be applied to this module as the next steps.

III. UAV DEMONSTRATION ARCHITECTURE

The Figure 1 shows our demonstration architecture of the cooperation between UAVs. This architecture consists of the Paparazzi modules in dark green. The paper focus on the Collision avoidance (using "TCAS") and Planner module that we have developed for our Paparazzi platform. "TCAS" is explained in the previous section and Planner is the mission planner module integrated in paparazzi.

Each aircraft is equipped with a Collision avoidance and an autopilot (AP) on board. The connection between GCS (Ground Control Station) and onboard aircraft is based on LINK module via Messages system and the IVY bus [12]. We send and bind messages on the IVY bus between Server, GCS and onboard aircraft. The data is sent with Messages are divided in three classes: Telemetry, Ground and Datalink. GCS sends Ground data to Aircraft. Server receives Telemetry raw data, dispatchs them and send them to GCS via Ground data. The mission will be defined by Mission definition module and then sent to Planner. The Planner interconnects to GCS for following aircrafts operations as interaction and monitoring level. Noting that IVY is an open-source software bus that can plug modules from many different languages such as C, C++, java, python, etc. Thus we can add and modify easily the new modules in paparazzi platform.

IV. SCENARIO

The first scenario we did for integrating the mission planner module in Paparazzi is shown in Figure 2. The scenario explains a common problem of configuration system in UAV domain. Multiple UAVs would cooperate together to solve separate tasks. The mission is composed of the surveillance task in a two-sections zone. For this mission, we have three UAVs (of the same type -fixed wing aircraft for the first simulation) and one operator. Basically, the main task is the surveillance of a large section with several UAVs. This task is decomposed to the surveillance task in two subsections each for only one UAV. The reconfiguration phase happens when an UAV is lacking of battery, the backup UAV will replace this one. The background tasks include then the separation phase -allocation of task and the deconflict phase -collision avoidance.

The goal of the mission planner is to allocate the two tasks to the three aircrafts. As the result (see Figure 3 (a) and (b) for more details), two tasks are allocated to two aircrafts. The other aircraft is circling outside around its standby point.

V. RESULTS

The simulation result is shown by the video https://www.dropbox.com/s/gy9djtkyrtqtdb7/video3AC1.avi in using Paparazzi simulation platform. Two surveillance sections are defined by S1-S2 and S2-S3 rectangle. We have three standby points for these three aircrafts. The HOME point is the point of starting our aircrafts and the other points are used for taking-off and landing. The mission planner will run as following steps:

• Launching three aircrafts: we launch our three aircrafts in blue, yellow and red (see Figure 4).

• Initializing planner: These three aircrafts are sent to their three standby points STDBY1, STDBY2, STDBY3. They are circling around these points and waiting for the next commands (see Figure 5).

• Starting planner: We are allocating two surveillance tasks to the three aircrafts (see Figure 6). Basically, they are distributed to the red and yellow aircraft because these two aircrafts are nearer from the two tasks than the blue aircraft. The trajectories of the surveillance task are defined by dividing the zone as the parallel lines. We will follow all these lines to cover all the zone. The radius of the circle once the aircraft goes out of the zone is determined in the flight plan before flying and can not be changed during flying (see Figure 7). At that moment, the status of two tasks is changed to ALLOCATED and the status of the red and yellow aircraft is changed to WORKING.

• Executing and replanning: Three aircrafts are continuing doing their current tasks. The planner is detecting any problem during flying including the battery, connection, etc. We simulate here a problem concerning the battery of the red aircraft by adding a "Kill AC" button. The red aircraft is supposed to land and we considered it to be dead (see Figure 8). The status of this aircraft is changed to DEAD. The planner detects the problem and make an online replanning. As the result, the blue aircraft will be sent to replace the red one (see Figure 9) (and the status of the blue aircraft is changed to WORKING after on). It is continuing the surveillance task of the red aircraft. We have now two working aircrafts are doing two tasks and one aircraft is dead (or preparing) outside. If we come to activate this aircraft again, we will send it to its standby point and waiting for the substitution.

On the computer Dell Latitude E6520 (Intel(R) Core(TM) i7-2640M CPU @ 2.80GHz, Ubuntu 12), the computation time measured during the simulation above are as following: We found that the computation time depends on the complexity of the mission such as the number of tasks, the number of aircrafts, the environment and the type of tasks, etc.

VI. CONCLUSION

The paper presents a new mission planner module integrated in Paparazzi platform. This planner allows us to allocate tasks to aircrafts in the context of cooperation between multiple UAVs. With this planner, we could carry out projects of surveillance tasks by Paparazzi. As the preliminary results of this paper, we perform a scenario of three aircrafts with the same type (fixed-wing aircraft). The simulation shows that we can deal with the replacement/substitution of multiple UAVs when there are problems of battery, connection missing, etc. The modality of the mission planner is currently based on the centralized communication system. For the next step of this work, we will simulate with other scenarios for different tasks, different types of aircrafts and use the ad-hoc network for the distributed communication system as well. We are going to apply this Mission planner paparazzi module in the D3CoS project to make a demonstrator of UAV domain. In the demonstrator, we will integrate not only our mission planner for the cooperation issue, but also the other tools to guarantee the human-machine interaction. The demonstration with real aircrafts will be performed at the end of this year.

  Task, Aicraft} Aircraft = {Aircraft Status, Sensor, Task} Task = {Task Position, Task Status, Goal} Waypoint = {x, y, altitude} Aircraft Status = {UNSTARTED, WORKING, ALIVE, DEAD} Task Status = {UNALLOCATED, ALLOCATED, DONE, FAILED}

Algorithm 1 :

 1 Task planning input : mission, task, aircraft data structure output: Allocating tasks to aircrafts // Loop for N tasks while j < N do if task[j].status = UNALLOCATED then // Plan for M UAVs for i ← 1 to M do if (aircraft[i].status = WORKING) & (aircraft[i].status = DEAD) then continue; end CostInsertTaskToPlan(aircraft[i],task[j]); MinimizeCost(aircrafts, task[j]); end InsertTaskToPlan(aircraft[i best],task[j]); end end Algorithm 2: Execution and task replanning input : mission, task, aircraft data structure output: Updating and Re-allocating tasks to aircrafts // Loop for N tasks while j < N do if task[j].status = DONE then // Replan for M UAVs for i ← 1 to M do if (aircraft[i].status = DEAD) & (task[j].status = FAILED) then RemoveTaskFromAircraft(aircraft[i], task[j]); end end end Do Algorithm 1 end Algorithm 3: Mission planner input : XML mission description data structure file output: Allocating tasks to aircrafts

Fig. 1 .

 1 Fig. 1. The architecture of our demonstration based on Paparazzi platform.

Fig. 2 .

 2 Fig. 2. The scenario with 3 aircrafts, 2 tasks and one operator.

Fig. 3 .

 3 Fig. 3. The mission planner process (a) and the distribution of two tasks to three aircrafts (b).

Fig. 4 .

 4 Fig. 4. Launching stage: three aircrafts in red, blue and yellow at the HOME point at the beginning. All the waypoints are in the UTM local coordinate relative to the the HOME point. We can see in the upper left corner the GCS interface of Paparazzi platform with strips of button for controlling/monitoring/editing before and during flight. The warning messages system is found on the lower left corner. The blocks of tasks and the altitude monitoring interface are on the lower side of the figure. At the end of the phase, three aircrafts are circling around the STDBY point.

Fig. 5 .

 5 Fig.5. Planner initialization: 3 aircrafts are circling around their standby points STDBY1, STDBY2, STDBY3 and waiting for the next commands. During flight, Collision avoidance module (using "TCAS") is always triggered each 500 (ms) and the conflict is detected if it happens. In this figure, we can see the conflict between the red and yellow aircraft. As the result, the red aircraft is climbing and the yellow one is descending for avoiding the collision between them. We can see this phase in the altitude interface on the lower side of the figure.

Fig. 6 .

 6 Fig.6. Start planner: the planner calculates and minimize the cost of the insertion of each task to each aircraft. Then, it chooses the yellow and the red aircrafts to two tasks of surveillance because they are nearest from the two tasks than the blue one.

Fig. 7 .

 7 Fig. 7. Start planner: two aircrafts are doing their surveillance tasks in two sections. We can see here the trajectory shapes during the surveillance task

Fig. 9 .

 9 Fig. 9. Updating planner: recalculating the cost of insertion and choosing the available blue aircraft to replace the red one for doing the S1-S2 section surveillance task

ACKNOWLEDGMENTS

This work is financed by the european project "D3CoS".