domain: math.math-ap nlin spi

text:
Andro Mikelić

Mary F Wheeler

Thomas Wick

Phase-field modeling of a fluid-driven fracture in a poroelastic medium *

Keywords: Hydraulic fracturing, Phase field formulation, Nonlinear elliptic-parabolic system, Computer simulations, Poroelasticity

In this paper we present a phase field model for a fluid-driven fracture in a poroelastic medium. In our previous work, the pressure was assumed given. Here we consider a fully coupled system where the pressure field is determined simultaneously with the displacement and the phase field. To the best of our knowledge, such a model is new in the literature. The mathematical model consists of a linear elasticity system with fading elastic moduli as the crack grows, which is coupled with an elliptic variational inequality for the phase field variable and with the pressure equation containing the phase field variable in its coefficients. The convex constraint of the variational inequality assures the irreversibility and entropy compatibility of the crack formation. The phase field variational inequality contains quadratic pressure and strain terms, with coefficients depending on the phase field unknown. We establish existence of a solution to the incremental problem through convergence of a finite dimensional approximation. Furthermore, we construct the corresponding Lyapunov functional that is linked to the free energy. Computational results are provided that demonstrate the effectiveness of this approach in treating fluid-driven fracture propagation.

Introduction

The coupling of flow and geomechanics in porous media is a major research topic in energy and environmental modeling. Of specific interest is induced hydraulic fracturing or hydrofracturing commonly known as fracking. This technique is used to release petroleum and natural gas that includes shale gas, tight gas, and coal seam gas for extraction. Here fracking creates fractures from a wellbore drilled into reservoir rock formations. According to the US National Petroleum Council 95% of oil and gas wells drilled today are hydraulically fractured. Clearly there are economic benefits of extracting vast amounts of formerly inaccessible hydrocarbons. In addition, there are environmental benefits of producing natural gas, much of which is produced in the United States from fracking. Opponents to fracking point to environmental impacts such as contamination of ground water, risks to air quality, migration of fracturing chemical and surface contamination from spills to name a few. For this reason, hydraulic fracturing is being heavily scrutinized resulting in the need for accurate and robust mathematical and computational models for treating fluid field fractures surrounded by a poroelastic medium.

Even in the most basic formulation, hydraulic fracturing is complicated to model since it involves the coupling of (i) mechanical deformation; (ii) the flow of fluids within the fracture and in the reservoir; (iii) fracture propagation. Generally, rock deformation is modeled using the theory of linear elasticity which can be represented by an integral equation that determines a relationship between fracture width and the fluid pressure. Fluid flow in the fracture is modeled using lubrication theory that relates fluid flow velocity, fracture width and the gradient of pressure. Fluid flow in the reservoir is modeled as a Darcy flow and the respective fluids are coupled through a leakage term. The criterion for fracture propagation is usually given by the conventional energy-release rate approach of linear elastic fracture mechanics (LEFM) theory; that is the fracture propagates if the stress intensity factor at the tip matches the rock toughness. Detailed discussion of the development of hydraulic fracturing models for use in petroleum engineering can be found in [START_REF] Adachi | Computer simulation of hydraulic fractures[END_REF] and in mechanical engineering and hydrology in [START_REF] Secchi | A method for 3-D hydraulic fracturing simulation[END_REF], [START_REF] De Borst | A numerical approach for arbitrary cracks in a fluidsaturated porous medium[END_REF], [START_REF] Irzal | A large deformation formulation for fluid flow in a progressively fracturing porous material[END_REF] and in references therein.

The goal of this paper is to develop a complete and robust phase field model for a fluid-driven fracture propagation in a poroelastic medium with the solid structure displacement, pressure and the phase-field variable as unknowns. We remark that in our preceding papers [START_REF] Mikelić | A phase-field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF], [START_REF] Mikelić | A quasistatic phase field approach to pressurized fractures[END_REF] and [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF], either pressure or displacements were assumed to be given.

The novelties in this coupled approach are:

• Fracture propagation in a poroelastic medium rather than in an elastic medium;

• Solving a fluid equation inside the fracture interacting with the fluid in the bulk poroelastic medium;

• Constructing the appropriate Lyapunov functional for the coupled three-field (displacementphase field-pressure) problem, which is linked to the phase-field modified free energy.

• Establishing mathematical stability of the three-field (displacement-phase field-pressure) equation, through constructing the corresponding Lyapunov functional and demonstrating effectiveness by providing several representative numerical tests.

The outline of our paper is as follows: In Section 2 we recall the basic equations of the sharp interface mesoscopic model and make some comments on the known phase-field approaches in the fracture mechanics. In Section 3, we present a detailed phase-field modeling of the quasistatic fracture propagation in a poroelastic medium. In Section 4 we introduce an incremental formulation of a phase-field model for a pressurized crack from which we derive a Lyapunov functional. In Section 5, we present a mathematical analysis of the incremental problem. In Section 6, a numerical formulation is briefly described. Finally in Section 7 we provide several numerical experiments demonstrating the capabilities of our proposed model.

Fundamental Equations

In order to fix ideas we address the simplest model of real applied importance, namely, the quasi-static single phase Biot system. We start by introducing the geometry.

Geometry

C Ω B = (0, L) 3 Figure 1: A crack C embedded in a porous medium. Here, the dimensions of the crack are assumed to be much larger than the pore scale size (black dots) of the porous medium.

Let C denote any open set homeomorphic to an ellipsoid in R 3 (a crack set). Its boundary is a closed surface ∂C. In most applications C is a curved 3D domain, with one dimension significantly smaller than the dominant two. An example is a penny shape crack.

The crack set C is surrounded by the poroelastic domain Ω = B \ C, where B = (0, L) 3 ⊃⊃ C. In many fracture propagation references, the crack C is considered to be a lower dimensional manifold and the lubrication theory is applied to describe the fluid flow (see e.g. [START_REF] Adachi | Computer simulation of hydraulic fractures[END_REF] and [START_REF] Ganis | Modeling fractures in a poro-elastic medium[END_REF]). We remark that 3D flow in C can be reconstructed from a lower dimensional lubrication approximation as described in [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF].

Nevertheless, we consider C as a 3D domain and use its particular geometry only when discussing the stress interface conditions.

The boundary of B = (0, L) 3

Biot equations

The quasi-static Biot equations (see e.g. [START_REF] Tolstoy | Acoustics, elasticity, and thermodynamics of porous media. Twenty-one papers[END_REF]) are an elliptic-parabolic system of PDEs, valid in the poroelastic domain Ω, where for every t ∈ (0, T ) we have

σ por -σ 0 = Ge(u) -αp R I; -div σ por = ρ b g; (1)
∂ t 1 M p R + div (αu) + div K R η R (ρ R g -∇p R ) = q, (2)
where σ 0 is the reference state total stress and q represents source terms. G is a constant symmetric positive definite rank-4 tensor. In our examples (7.1) and ( 7.2) we model the stress Ge(u) using a linear stress-strain relationship:

Ge(u) = 2µe(u) + λtr(e(u))I,
where I is the identity matrix, tr( The important parameters and unknowns are listed in the Table 1.

Due to the fracture size, the relative velocity satisfies Darcy's law

v F = K F η F (G -∇p F ). (3)
Note that for a penny shaped C the dominant velocity part is two-dimensional. An explicit expression for the permeability K F as a function of the fracture shape is derived in [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF]. The term G is a linear function of the gravity and leakage (see [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF], Sec 2.3 for details). In the fracture C, the pressure satisfies the mass conservation equation

c F ∂ t p + div K F η F (G -∇p) = q. ( 4
)
Modeling the moving interface ∂C is delicate. We impose the following interface conditions p R = p F and σ por n = -p F n on the moving interface ∂C (5)

v F • n = v R • n = K R η R (ρ R g -∇p R ) • n on the moving interface ∂C, ( 6
)
where n is the normal exterior to the fracture. Condition ( 6) is discussed in [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF], Sec 2.3. Condition ( 5) is a simplification of the pressure interface conditions and more careful modeling by upscaling indicates presence of an additional term. We do not dwell on this question here.

We are primarily interested in a fracture propagating due to pressure (water) injection effects. Consequently, we keep the source/sink terms q and set σ 0 = 0 and neglect gravity effects.

Remark 1. The main difference in the present work as compared to [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF] is that in [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF], we solved the (u, ϕ) system monolithically and coupled this system through fixed-stress to the pressure diffraction equation in a sequential procedure. In the present paper, we formulate a fully-coupled system that now allows us to formulate a free energy functional.

The strategy developed in [START_REF] Mikelić | A quasistatic phase field approach to pressurized fractures[END_REF]- [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF] relies on regularity. For calculations of the displacement and phase field in [START_REF] Mikelić | A quasistatic phase field approach to pressurized fractures[END_REF], we required p ∈ W 1,1 (0, T ; W 1,2r ((0, L) 3 )), with r > 3. In general, it is not clear that such smoothness can be guaranteed. In [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF] we assumed a sharp boundary between the poroelastic domain and the fracture for which C 1 regularity was needed. The new approach here does not require higher regularity of the phase field. Hence it can be used as a general formulation. The strategy applied in [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF] can then be applied to obtain a more accurate approximation of the interface between the poroelastic medium and the fracture.

Modeling fractures with a phase-field approach in solid mechanics

In the following, we present a continuous three-field hydraulic fracturing model. Then we discretize in time and provide a complete mathematical theory for the equations of the corresponding incremental formulation of the system. The mathematical model involves the coupling of a linear elasticity system with a) the pressure equation and b) with an elliptic variational inequality for the phase field variable. With this approach, branching of fractures and heterogeneities in mechanical properties can be effectively treated as demonstrated numerically in Section 7.

For the elasticity part, our formulation follows Miehe et al. in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[END_REF] and enters a generalization of the thermodynamically consistent framework for phase-field models of quasi-static crack propagation in elastic solids, together with incremental variational principles. The phase field equations are based on the variational approach to elastic fractures by Francfort and Marigo ([9] and [START_REF] Bourdin | The Variational Approach to Fracture[END_REF]). Our methodology represents an extension to cracks in a poroelastic medium containing a viscous fluid.

In fracture mechanics, Griffith's criterion states that crack propagation occurs when the elastic energy restitution rate reaches its critical value G c . If τ is the traction force applied at the part of the boundary ∂ N Ω, then one can associate to the crack C the following total energy

E(u, C) = Ω 1 2
Ge(u) : e(u) dx -

∂ N Ω τ • u dS - Ω F div u dx + G c H 2 (∂C), (7)
where F is a volume force. In fracture mechanics the crack is considered as a lower dimensional manifold and would correspond to ∂C.

In [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and [START_REF] Bourdin | The Variational Approach to Fracture[END_REF] a phase field regularization of the term G c H 2 (∂C) was used. The equations in the displacement unknown u and the phase-field unknown ϕ are derived as the corresponding Euler-Lagrange equations. They are used for the determination of the cracks spreading in fracture mechanics.

The situation in poroelasticity is more complicated than the classical one from solid mechanics.

Specific challenges of phase-field fracture in poroelasticity

First, we recall that our fracture is a (thin) 3D body and not a surface, with dimensions being small compared to B, but large compared to the characteristic pore size . At the pore level we have a complex interaction between fracture and the pore skeleton, but such a modeling is presently out of reach. Nevertheless, the experimental observation is that if the stresses (respectively the free energy) cross a threshold, the fracture will start growing and propagating. We will follow the ideas from [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and [START_REF] Bourdin | The Variational Approach to Fracture[END_REF] and add to the free energy the phase field regularization of the term G c H 2 (∂C). Since the observed and simulated fractures correspond to smeared two-dimensional surfaces, our treatment is consistent with ideas from the literature (see e.g. [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[END_REF]).

Since we do not have a total energy functional for Biot's system (1)-( 2), the energy functional (7) must be modified to include pressure effects and a term describing the propagation for some threshold value of the stress. Our strategy is to generalize the notion of the free energy for poroelasticity in order to include the phase-field. This is achieved by regularizing the sharp fracture surface topology in the poroelastic medium by diffusive crack zones described by a scalar auxiliary variable. This variable is a phase-field that interpolates between the unbroken and the broken states of the material. After making all these modifications and taking into the account the pressure field, we will not have an energy functional; however we will derive a generalized free energy. The latter will allow us to construct a Lyapunov functional and as a consequence the incremental problem is well-posed.

3 Formulation of a fully-coupled displacement, phase-field, pressure system

In fracture mechanics the phase field unknown ϕ is defined on B × (0, T ) by the regularized crack functional

Γ ε (ϕ) = G c B ( 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 ) dx = G c B γ(ϕ, ∇ϕ) dx, ( 8
)
where γ is the crack surface density per unit volume and ε a positive constant. This regularization of H 2 (∂C), in the sense of the Γ-limit when ε → 0, was used in [START_REF] Bourdin | The Variational Approach to Fracture[END_REF]. Even if we do not deal with a surface and its energy, it is observed that the fracture grows in a porous medium if the stress is greater than some critical threshold G c . Here, G c depends on the skeleton compressibility and on the fluid viscosity. In classical fracture mechanics setting G c corresponds to the energy release rate. In deference to classical fracture mechanics, the meaning of the problem for ε = 0 is not clear and our ε will be a small positive parameter, but we will not aim to pass to the limit ε → 0. The crack phase-field unknown ϕ is introduced through a regularized elastic energy functional minimization for given pressure field. Then the form of Γ ε (ϕ) yields that meas{(x, t) ∈ B ×(0, T ) | ϕ(x, t) ≤ a < 1} ≤ Cε. Therefore the parameter ε represents the characteristic length scale and formally the crack behaves as a surface when ε → 0.

Our further considerations are based on the fact that evolution of cracks is fully dissipative in nature. First, the crack phase field ϕ is intuitively a regularization of 1 -χ C and we impose its negative evolution

∂ t ϕ ≤ 0. (9)
After the seminal work by Biot [START_REF] Biot | Mechanics of Deformation and Acoustic Propagation in Porous Media[END_REF], the Biot stress and the pressure are determined using the strain energy W of a porous elastic medium. W is the free isothermal energy density of the fluid-solid system. After [START_REF] Biot | Mechanics of Deformation and Acoustic Propagation in Porous Media[END_REF], W depends on the strain tensor e(u) and on the increment of the fluid content

ζ = p/M + α div u, i.e. W = W (e, ζ). Since it is required that p = ∂W ∂ζ and σ = ∂W ∂e ,
we have

W (e, ζ) = 1 2 Ge(u) : e(u) + p 2 /(2M ) = 1 2 Ge(u) : e(u) + M 2 ζ -αdiv u 2 .
Our idea is to generalize this energy and to use the 'phase field' free energy. It has to take into the account the following effects:

• The degradation of the elastic properties in the fracture, described by the multiplication of the Gassmann tensor G by the function A, such that A(ϕ) is a smooth monotone increasing function satisfying

A(0) = k > 0, k ε, A (0) = 0 and A(1) = 1. ( 10
) A usual choice is A(ϕ) = (1 -k)ϕ 2 + k.
• The change of the porosity being different in the poroelastic part and in the fracture, and now described through a function c satisfying

c(ϕ) is a smooth bounded function such that c(0) = c F > 0, inf z c(z) = c min > 0, c (0) = 0 and c(1) = 1/M. ( 11
)
Note that the coefficients A and c are naturally defined on [0, 1]. We extend them for negative values of argument by the value at ϕ = 0 and pose A(ϕ) = A(ϕ + ) and c(ϕ) = c(ϕ + ). Due to constraint (9), ϕ will always remain smaller or equal to 1. For more discussion we refer to [START_REF] Mikelić | A quasistatic phase field approach to pressurized fractures[END_REF] and further comments bellow. Our new 'phase-field' free energy density including fractures is now

W ε (e, ζ, ϕ) = 1 2 A(ϕ + )Ge(u) : e(u) + c(ϕ + )p 2 2 + G c 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 .
We note the presence of the term allowing the fracture to grow, once the threshold energy is crossed. Furthermore the phase field unknown ϕ should allow a smooth change from the fracture zone to the poroelastic zone. For the fracture spreading it is crucial to accurately model the interface and it has to enter the phase field equation on the correct way.

Our strategy is to generalize the well-established approach from fracture mechanics. We proceed in several steps:

The elastic part energy functional

We suppose a given pressure, but bear in mind that this pressure is linked to the displacement and the phase-field unknown. Then we borrow the energy functional from fracture theory of solid mechanics and modify by taking into the account the porosity change. Following modeling of a brittle fracture of fluid-infiltrated material from [7], Sec. 4.4.2., we add a cross term involving the integral of the product between the pressure and the porosity change. This is the main difference with the approaches in the simplified case [START_REF] Mikelić | A phase-field approach to the fluid filled fracture surrounded by a poroelastic medium[END_REF] , [START_REF] Mikelić | A quasistatic phase field approach to pressurized fractures[END_REF] and [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF]. This yields the following functional:

Ẽε (u, ϕ) = B 1 2 A(ϕ)Ge(u) : e(u) dx - ∂ N B τ • u dS - ∂C τ • u dS - B αϕpdiv u dx + G c B 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 dx + B 1 2 ( 1 M -c(ϕ + ))p 2 dx. (12)
The last term is constant with respect to the displacement u and does not contribute to the displacement equation.

Interface law for the normal stresses on the fracture boundary

Now we observe that the interface term ∂C τ • u dS (12) has to be transformed into a volume term.

Here, we model the fracture pressure using an interface law. A general description of a crack embedded in a porous medium is illustrated in Figure 1. We concentrate on a simplified situation in which the complex interaction crack/pore structure is neglected. We recall that the crack is a 3D thin domain with width much smaller than length, then lubrication theory can be applied. Hence, the leading order of the stress in C is -p f I, where p f denotes the fracture pressure. From the poroelastic side the stress is given by σ = σ 0 + Ge(u) -αp B I, where p B denotes the Biot pressure in the poroelastic medium. At the crack boundary we have the continuity of the contact force, which yields

σn = (Ge(u) -αp B I)n = -p f n, ( 13
)
where n is the normal vector exterior to the poroelastic domain.

Next, we eliminate the traction crack surface integrals and obtain

∂C σnw dS = - ∂C p f w n dS = - Ω div (p B w) dx + ∂ N B p B w n dS = - Ω p B div w dx - Ω ∇p B w dx + ∂ N B p B w n dS,
where w n denotes the normal component of the vector function w.

Remark 2. In the previous calculations, the first line reads

∂C σnw dS = - ∂C p f w n,poroelastic dS.
Setting now p = p f in the crack and p = p B in the poroelastic medium, then the above calculations yield

- ∂C p f w n,poroelastic dS = C ∇pw dx - C div (pw) dx - ∂C p f w n,poroelastic dS = C ∇pw dx - ∂C p(w n,poroelastic -w n,crack ) dS.
The last term coincides with the virtual work of the pressure force as introduced in [START_REF] Calhoun | A two dimensional asymmetrical crack problem[END_REF].

In the above calculations, surface integrals are now treated with Gauss' divergence theorem:

-

∂ N Ω τ • w dS + ∂C pw n dS - Ω αp div w dx = - Ω (α -1)p div w dx+ Ω ∇pw dx - ∂ N Ω τ • w dS = - Ω (α -1)p div w dx+ Ω ∇p • w dx - ∂ N Ω (τ • w + pw n ) dS. (14)
Using [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF], Ẽε transforms to

Êε (u, ϕ) = B 1 2 A(ϕ + )Ge(u) : e(u) dx - ∂ N B (τ • u + pu n ) dS + Ω ((1 -α)pdiv u + ∇p • u) dx +G c B 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 dx + B ( 1 M - c(ϕ + ) 2 )p 2 dx. ( 15
)
We suppose that the crack C does not interact with ∂ N Ω. Hence there is no need to insert ϕ in the boundary integral over ∂ N B.

Phase-field formulation of the pressure cross terms

It remains to extend the pressure cross term

Ω ((1 -α)pdiv u + ∇p • u) dx
to a term with integration over B by adding ϕ. Due to the complexity of our model, it is not clear that we will be able to prove the nonnegativity of ϕ. Changing the pressure cross term to

B ϕ (1 -α)pdiv u + ∇p • u)
dx is likely to give nonzero values in the domain where ϕ is negative. This is not admissible. Using ϕ + instead of ϕ is also unsatisfactory, since the energy functional would not be C 1 any more. The lack of regularity would lead to complications in numerical simulations. Our approach is to use the function ϕ 2 + instead of ϕ. For 0 ≤ ϕ ≤ 1, the choice between ϕ 2 + and ϕ + in the pressure cross term does not affect the phase field approximation. If ϕ = 1 -χ C , the two functions coincide.

After the above transformations, we have the following variant of the energy functional ( 15):

E ε (u, ϕ) = B 1 2 A(ϕ + )Ge(u) : e(u) dx - ∂ N B (τ • u + pu n ) dS + B ϕ 2 + (1 -α)pdiv u + ∇p • u) dx +G c B 1 2ε (1 -ϕ) 2 + ε 2 |∇ϕ| 2 dx + B ( 1 M - c(ϕ + ) 2 )p 2 dx. (16)

Fréchet derivative of E ε

As in the literature, {u, ϕ} are critical points of the non-convex functional E ε . Taking the Fréchet derivatives yields

-div A(ϕ + )Ge(u) + ϕ 2 + ∇p + (α -1)∇(ϕ 2 + p) = 0 in B, (17)
L(u, ϕ) = -G c ε∆ϕ - G c ε (1 -ϕ) + A (ϕ + ) 2 Ge(u) : e(u) - c (ϕ + ) 2 p 2 + 2(1 -α)ϕ + p div u + 2ϕ + ∇p • u ≤ 0 in B, (18)
∂ t ϕ ≤ 0 in B, (19)
∂ t ϕ • L(u, ϕ) = 0 in B. (20)
Here, the last expression ( 20) is a complementarity condition (namely the Rice condition) linked to the irreversibility constraint [START_REF] Schrefler | On adaptive refinement techniques in multi-field problems including cohesive fracture[END_REF].

Adding reservoir and fracture pressure equations

We have seen that under the assumption of the weak compressibility of the fluid, the mass conservation equation reads [START_REF] Bangerth | Differential Equations Analysis Library[END_REF] in the poroelastic part and (4) in the fracture. Following [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF] it can be written as the pressure diffraction problem

θ∂ t p -div K ef f (∇p -G) = q in B, (21)
where

θ = θ(x, t) = χ Ω 1 M + χ C c F , q = q(x, t) := χ Ω q R ρ R -αχ Ω ∂ t div u + χ C q F ρ F ; G = χ C G, K ef f = χ Ω K R η R + χ C K F η F .
We recall that gravity effects are neglected.

The formulation contains the interface between Ω and C, which changes in time, as the crack propagates. Our idea is to replace the diffraction problem ( 21) by a phase field equation for the pressure diffraction equation. Following the distinction between terms containing α -1 as factor and others from equation ( 17), we have the following global problem for the pressure:

∂ t (c(ϕ + )p) -div K ef f (ϕ + )(∇p -G(ϕ + )) + (α -1)∂ t (ϕ 2 + div u)+ div ∂ t (ϕ 2 + u) = q in B. (22)
Note that for ϕ = 0 and ϕ = 1 we have again equation [START_REF] Tolstoy | Acoustics, elasticity, and thermodynamics of porous media. Twenty-one papers[END_REF]. G(0) = G and G(1) = 0.

The fully coupled system for {u, ϕ, p}

We now summarize the Euler-Lagrange equations from STEP 4 and the phase-field pressure diffraction problem from STEP 5 and obtain the following fully-coupled system: Find u, ϕ, p such that:

-div A(ϕ + )Ge(u) + ϕ 2 + ∇p + (α -1)∇(ϕ 2 + p) = 0 in B, (23)
u = 0 and -K ef f (ϕ + )(∇p -G(ϕ + )) • n = v inj on ∂ D B, (24)
-Ge(u

)n = τ + p bdry n on ∂ N B, (25)
-G c ε∆ϕ - G c ε (1 -ϕ) + 1 2 A (ϕ + )Ge(u) : e(u) - 1 2 c (ϕ + )p 2 +2ϕ + ((1 -α)p div u + ∇p • u ≤ 0 in B, (26)
∂ t ϕ ≤ 0 on B (27) -G c ε∆ϕ - G c ε (1 -ϕ) + 1 2 A (ϕ + )Ge(u) : e(u) - 1 2 c (ϕ + )p 2 +2ϕ + ((1 -α)p div u + ∇p • u ∂ t ϕ = 0 in B, (28)
ϕ(x, 0) = ϕ 0 (x) and ∂ϕ ∂n = 0 on ∂B;

(29)

∂ t (c(ϕ + )p) -div K ef f (ϕ + )(∇p -G(ϕ + )) + (α -1)∂ t (ϕ 2 + div u) +div ∂ t (ϕ 2 + u) = q in B, (30) p
(x, 0) = p 0 (x) on B; p = p bdry on ∂ N B, (31)
where ( 24) and (31) contain the boundary conditions for the displacement and the pressure.

Weak formulation of the fully coupled system

As functional spaces of admissible displacements and pressures, we choose

V U = {z ∈ H 1 (B) 3 | z = 0 on ∂ D B } V P = {π ∈ H 1 (B) | π = 0 on ∂ N B }.
We proceed by writing the variational formulation corresponding to problem ( 23)-(31): Find {u, ϕ, p} such that

B A(ϕ + )Ge(u) : e(w) dx + B ϕ 2 + ∇p • w dx + (α -1) B ∇(ϕ 2 + p) • w dx = ∂ N B (τ • w + p bdry w n ) dS, ∀w ∈ V U a.e. on (0, T ); (32) B(ϕ, ψ) = G c B ε∇ϕ • ∇ψ - 1 ε (1 -ϕ)ψ dx + 1 2 B A (ϕ + )Ge(u) : e(u) -c (ϕ + )p 2 ψ dx + 2 B ϕ + ((1 -α)p div u + ∇p • u ψ dx ≤ 0 a.e. on (0, T ), ∀ψ ∈ H 1 (B), ψ ≥ 0 a.e. on B; (33)
∂ t ϕ ≤ 0 a.e. on B × (0, T );

(34)

B(ϕ, ∂ t ϕ) = 0 a.e. on (0, T ); (35) ϕ(x, 0) = ϕ 0 (x) on B, 0 ≤ ϕ 0 (x) ≤ 1; (36) d dt B c(ϕ + )pz dx + B K ef f (ϕ + )(∇p -G(ϕ + )) • ∇z dx+ (α -1) B ∂ t (ϕ 2 + div u)z dx + B div (∂ t (ϕ 2 + u))z dx = B qz dx - ∂ D B
v inj z dS, ∀z ∈ V P a.e. on (0, T ); (37)

p(x, 0) = p 0 (x) on B; p = p bdry on ∂ N B. (38)
We notice that B(ϕ, ∂ t ϕ) = 0 is just another notation for the complementarity condition.

We will undertake the study of the incremental (i.e. the time discretized problem) in Sec. 4. Nevertheless, the basic idea of this phase-field formulation is to formulate a Lyapunov functional for our full model that generalizes the free energy. The corresponding calculations follow in the next step.

Construction of the Lyapunov functional

+ G c ε|∇ϕ| 2 + 1 ε (1 -ϕ) 2 + c(ϕ + )p 2 - 2c(ϕ + )pp bdry -2(α -1)ϕ 2 + div u p bdry -2 div (ϕ 2 + u) p bdry dx -2 ∂ N B τ • u + (1 -α)p bdry u n dS + B K ef f (ϕ + )|∇p| 2 dx = - ∂ N B ∂ t τ • u + ∂ t p bdry u n dS - B K ef f (ϕ + )(∇p -G(ϕ + )) • ∇p bdry dx - B c(ϕ + )p∂ t p bdry dx -(α -1) B ϕ 2 + div u ∂ t p bdry dx - B div (ϕ 2 + u) ∂ t p bdry dx + B K ef f (ϕ + ) G(ϕ + ) • ∇p bdry dx - ∂ D B v inj (p -p bdry ) dS + B q(p -p bdry ) dx. ( 39
)
Proof. We will discuss the proof only for the incremental problem, in Theorem 2.

Remark 3. The free energy for problem (32)-( 38) is

{u, ϕ, p} → 1 2 A(ϕ + )Ge(u) : e(u) + G c ε|∇ϕ| 2 + 1 ε (1 -ϕ) 2 + c(ϕ + )p 2 . ( 40
)
4 An incremental formulation of the fully-coupled system

In the following, we consider a quasi-static formulation where velocity changes are small. First, we derive an incremental form, i.e., we replace the time derivative in inequality ( 9) with a discretized version; more precisely

∂ t ϕ → ∂ ∆t ϕ = (ϕ -Φ)/(∆t),
where ∆t > 0 is the time step and Φ is the phase field from the previous time step. After time discretization, our quasistatic constrained minimization problem becomes a stationary problem, called the incremental problem.

The entropy condition ( 9) is imposed in its discretized form and we introduce a convex set K:

K = {ψ ∈ H 1 (B) | ψ ≤ Φ ≤ 1 a.e. on B}. ( 41
)
Note that the value of the phase field unknown ϕ from the previous time step enters only the convex set K, as the obstacle Φ. U and P will denote the values of the displacement and the pressure from the previous time step.

Let us write the incremental form of problem ( 23)-( 31).

-div A(ϕ + )Ge(u) + ϕ 2 + ∇p + (α -1)∇(ϕ 2 + p) = 0 in B, (42)
u = 0 and -K ef f (1)(∇p -G(Φ + )) • n = v inj on ∂ D B, (43)
Ge(u)n = τ + p bdry n on ∂ N B, (44)
-G c ε∆ϕ - G c ε (1 -ϕ) + 1 2 A(ϕ)Ge(U) : e(U) - 1 2 J (ϕ)P (2p -P ) +D(ϕ) (1 -α)p div U + ∇p • U ≤ 0 in B, (45)
∂ ∆t ϕ ≤ 0 on B (46) -G c ε∆ϕ - G c ε (1 -ϕ) + 1 2 A(ϕ)Ge(U) : e(U) + 1 2 J (ϕ)P (2p -P ) +D(ϕ) ((1 -α)p div U + ∇p • U ∂ ∆t ϕ = 0 in B, (47)
ϕ(x, 0) = ϕ 0 (x) and ∂ϕ ∂n = 0 on ∂B;

(48)

∂ ∆t (c(ϕ + )p) -div K ef f (Φ + )(∇p -G(Φ + )) + (α -1)∂ ∆t (ϕ 2 + div u) +div ∂ ∆t (ϕ 2 + u) = q in B, (49)
ϕ(x, 0) = ϕ 0 (x) and p(x, 0) = p 0 (x) on B; p = p bdry on ∂ N B,

where (47) is the strong form of Rice' condition. The functions A, J and D correspond to discrete derivatives of A, c and ϕ 2 and are given by

A(ϕ) = A(ϕ + ) -A(Φ + ) ϕ -Φ ; J (ϕ) = c(ϕ + ) -c(Φ + ) ϕ -Φ and D(ϕ) = ϕ 2 + -Φ 2 + ϕ -Φ . ( 51
)
Here, the system (42)-( 47) is a variational inequality; in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[END_REF] a penalization term is used for solving the inequality. To the best of our knowledge the modeling of the full system and its analysis are new in the literature.

Well-posedness of the model

We start by making the following assumptions on the data Hypothesis 1.

G is a positive definite constant rank-4 tensor, G is a bounded continuous function of ϕ + , x and t, K ef f is a continuous bounded symmetric matrix, positively definite, with the ellipticity constant independent of the argument, A, c are continuously differentiable functions, taking values between two positive constants;

Φ ∈ H 1 (B), Φ ≤ 1 a.e. on B; τ ∈ L 2 (∂ N B)), p bdry ∈ H 1 (B), v inj ∈ L 2 (∂ D B) and q ∈ L 2 (B).

A finite dimensional approximation

Let {ψ r } r∈N be a smooth basis for H 1 (B), {π r } r∈N a smooth basis for V P and {w r } r∈N be a smooth basis for V U . We start by defining a finite dimensional approximation to the problem (42)-(50). ..,M = {ξ1 , ξ 2 , ξ 3 } and X the finite dimensional space spanned by the set of all such ξ. X is isomorphic to R 3M and we take the natural scalar product. After setting

P 1,r (ξ) = B A( φM )Ge(u M ) : e(w r ) + (α -1)∇(( φM ) 2 p M ) • w r dx+ B ( φM ) 2 ∇p M • w r dx - ∂ N B (τ • w r + p bdry w r n ) dS, r = 1, . . . , M ; (57)
P 2,r (ξ) = G c B - 1 ε (1 -ϕ M )ψ r + ε∇ϕ M • ∇ψ r dx+ B ψ r δ(ϕ M -Φ) + dx + B D( φM ) ((1 -α)p M div U + ∇p M • U ψ r dx + 1 2 B A( φM )Ge(U) : e(U) -J ( φM )P (2p M -P ) ψ r dx, ∀ r = 1, . . . , M, (58)
P 3,r (ξ) = B c( φM )p M + div (( φM ) 2 u M ) + (α -1)( φM ) 2 div u M π r dx +∆t B K ef f (Φ + )(∇p M -G(Φ + )) • ∇π r dx -∆t B qπ r dx+ ∆t ∂ D B v inj π r dS + B c(Φ + )P + div (Φ 2 + U)+ (α -1)Φ 2 + div U π r dx, ∀ r = 1, . . . , M, (59)
we see that problem (52)-(54) has a solution if and only if equation P(ξ) = 0 has a solution.

The nonlinear mapping P is obviously continuous between X and X. Using a well-known corollary 1 of Brouwer's fixed point theorem, it is enough to prove that (P(ξ), ξ) X > 0 for ξ with sufficiently large norm. Existence of at least one root would follow.

We start by multiplying P 1,r (ξ) by a r and taking the sum with respect to r. It yields

(P 1 (ξ), ξ 1 ) = B A( φM )Ge(u M ) : e(u M ) + (α -1)∇(( φM ) 2 p M ) • u M dx + B ( φM ) 2 ∇p M • u M dx - ∂ N B (τ • u M + p bdry u M n ) dS. (60)
Since u M ∈ V U , the trace theorem and Korn's inequality yield

| ∂ N B (τ • u M + p bdry u M n ) dS| ≤ C||e(u M )|| L 2 (B) (||τ || L 2 (∂ N B) 3 +||p bdry || L 2 (∂ N B) ); (61) A( φM )Ge(u M ) : e(u M ) ≥ C K ||u M || H 1 (B) . (62)
and in (P 3 (ξ), ξ 3 ). Putting them together yields

B (α -1) ∇(( φM ) 2 (p M -p bdry )) • u M + ( φM ) 2 div u M (p M -p bdry ) dx = (α -1) B div ( φM ) 2 (p M -p bdry )u M dx = 0; B ( φM ) 2 ∇(p M -p bdry ) • u M + div (( φM ) 2 u M )(p M -p bdry ) dx = 0.
Therefore we have

(P(ξ), ξ) = (P 1 (ξ), ξ 1 ) + (P 2 (ξ), ξ 2 ) + P 3 (ξ), ξ 3 ) ≥ G c B (ϕ M ) 2 4ε + ε|∇ϕ M | 2 dx +C||u M || 2 H 1 (B) 3 + C||p M || 2 H 1 (B) -C||P || 2 H 1 (B) -C||U|| 2 H 1 (B) 3 -C -C||Φ|| H 1 (B) 3 + δ B (ϕ M -Φ) 2 + dx. (68)
It follows from ( 68) that (P(ξ), ξ) > 0 for ||ξ|| = R, with sufficiently large R.

Obviously corresponding solutions {u M , ϕ M , p M } satisfy a priori estimate (56).

Existence of a solution to the incremental problem

In this section we show existence of a weak solution:

Theorem 1. Assume the Hypothesis 1. Then there exists at least one variational solution {u, ϕ, p} ∈ V U × H 1 (B) ∩ K × V P for problem (42)-(50).

Proof. (of Theorem 1) By Proposition 2 there is a solution {u M , ϕ M , p M } for problem (52)-(54), satisfying a priori estimate (56). Therefore there exists {u, ϕ, p} and a subsequence, denoted by the same superscript, such that

{u M , ϕ M , p M -p bdry } → {u, ϕ, p -p bdry } weakly in V U × H 1 (B) × V P ,
strongly in L q (B) 5 , q < 6, and a.e. on B, as M → ∞.

Obviously (ϕ M -Φ) + → 0, as M → ∞, and ϕ ∈ K.

Since φM → φ strongly in L q B, for all q < +∞, we pass to the limit in equation (52) without obstacles. Therefore, the triple {u, ϕ, p} satisfies equation (42) and boundary condition (44).

The analogous conclusion holds for variational equation (54) and the triple {u, ϕ, p} satisfies equation (49) and boundary conditions (43), (50).

Next, let ψ ∈ K N + = {z ∈ W N : z(x) ≤ 0 a.e. on B}, N ≤ M . Then we have

G c B - 1 ε (1 -ϕ M )ψ + ε∇ϕ M • ∇ψ dx + B M (ϕ M -Φ) + ψ dx + B D( φM ) (1 -α)p M div U + ∇p M • U ψ dx+ 1 2 B A( φM )Ge(U) : e(U) -J ( φM )P (2p M -P ) ψ dx = 0. ( 70
)
Passing to the limit in equation ( 70) is straightforward and we conclude that

G c B - 1 ε (1 -ϕ)ψ + ε∇ϕ • ∇ψ dx + B D( φ) (1 -α)p div U + ∇p • U ψ dx+ 1 2 B A( φ)Ge(U) : e(U) -J ( φ)P (2p -P ) ψ dx ≥ 0. ( 71
)
After passing to the limit N → ∞, we conclude that {u, ϕ, p -p bdry } ∈ V U × H 1 (B) ∩ K × V P is a solution to the inequality (45) and ( 46) is satisfied. It remains to prove the Rice condition (47).

To prove it, we observe that with the choice ψ = ϕ-Φ inequality (71) holds true. Finally, we choose ϕ M -Φ as a test function in equation (70). Using the weak lower semi-continuity of the quadratic terms yields

G c B - 1 ε (1 -ϕ)(ϕ -Φ) + ε∇ϕ • ∇(ϕ -Φ) dx + B D( φ) (1 -α)p div U + ∇p • U (ϕ -Φ) dx+ 1 2 B A( φ)Ge(U) : e(U) -J ( φ)P (2p -P ) (ϕ -Φ) dx ≤ 0. (72)
After comparing (71) and (72), we obtain Rice's condition (47).

The Lyapunov functional for the incremental problem

In this subsection we construct the Lyapunov functional corresponding to the sequence of incremental problems (42)-(50).

We suppose the quasistatic problem ( 23)-( 31) is discretized with a uniform time step ∆t. Given solutions at discrete times t j , j = 0, . . . , N , {u j , ϕ j , p j }, with Φ = ϕ j-1 , U = u j-1 and P = p j-1 , are extended from the discrete times {t j } 0≤j≤N to (0, T ) by

ϕ ∆t (t) = ϕ j if t j ≤ t < t j+1 , j = 0, . . . N -1;
(73)

ϕ ∆t (t 0 ) = ϕ 0 ; (74) u ∆t (t) = u j if t j ≤ t < t j+1 , j = 0, . . . N -1. ( 75
)
p ∆t (t) = p j if t j ≤ t < t j+1 , j = 0, . . . N -1; (76) p ∆t (t 0 ) = p 0 . ( 77
)
u ∆t (t 0 ) ∈ V U is the solution for (42), ( 43), (44). In the formulation we need values at t -∆t, which correspond to the previous time step. We set

Φ ∆t (t) = ϕ ∆t (t -∆t); U ∆t (t) = u ∆t (t -∆t) and P ∆t (t) = p ∆t (t -∆t) ∆t ≤ t ≤ T.
The variational formulation of the incremental problem is then Find {u ∆t , ϕ ∆t , p ∆t } ∈ V U × H 1 (B) × V P such that ϕ ∆t (t) ≤ ϕ ∆t (t -∆t) for every t ∈ (0, T ) and

B A( φ∆t )Ge(u ∆t ) : e(w) + ( φ∆t ) 2 ∇p ∆t + (α -1)∇(( φ∆t ) 2 p ∆t ) • w dx = ∂ N B (τ ∆t • w + p bdry ∆t w n ) dS, ∀w ∈ V U , (78)
G c B - 1 ε (1 -ϕ ∆t )ψ + ε∇ϕ ∆t • ∇ψ dx + B D( φ∆t ) ((1 -α)p ∆t div U ∆t + ∇p ∆t • U ∆t ψ dx+ 1 2 B A( φ∆t )Ge(U ∆t ) : e(U ∆t ) -J ( φ∆t )P ∆t (2p ∆t -P ∆t ) ψ dx ≤ 0, ∀ψ ∈ H 1 (B), ψ ≥ 0, a.e. on B, (79)
G c B - 1 ε (1 -ϕ ∆t )(ϕ ∆t -Φ ∆t ) + ε∇ϕ ∆t • ∇(ϕ ∆t -Φ ∆t ) dx + B D( φ∆t ) ((1 -α)p ∆t div U ∆t + ∇p ∆t • U ∆t (ϕ ∆t -Φ ∆t ) dx+ 1 2 B A( φ∆t )Ge(U ∆t ) : e(U ∆t ) -J ( φ∆t )P ∆t (2p ∆t -P ∆t ) (ϕ ∆t -Φ ∆t ) dx = 0, (80)
B ∂ ∆t c( φ∆t )p ∆t + div (( φ∆t ) 2 u ∆t ) + (α -1)( φ∆t ) 2 div u ∆t π dx

+ B K ef f ( Φ∆t ) ∇p ∆t -G(Φ ∆t ) • ∇π dx = B q ∆t π dx - ∂ D B v inj ∆t π dS, ∀ π ∈ V P , (81)
with

A( φ∆t ) = A( φ∆t ) -A( Φ∆t ) ϕ ∆t -Φ ∆t ; J ( φ∆t ) = c( φ∆t ) -c( Φ∆t ) ϕ ∆t -Φ ∆t and D( φ∆t ) = φ2 ∆t -Φ2 ∆t ϕ ∆t -Φ ∆t . ( 82
)
We start with several auxiliary results, which will allow us to construct the Lyapunov functional:

Lemma 1. Let N ∈ N. Then we have ∆t N -1 n=0 A( φn+1 )Ge(u n+1 ) : e(∂ ∆t u n+1 ) = 1 2 A( φN )Ge(u N ) : e(u N )- A( φ0 )Ge(u 0 ) : e(u 0 ) + 1 2 N -1 n=0 A( φn+1 )Ge(u n+1 -u n ) : e(u n+1 -u n )- ∆t 2 N -1 n=0 A( φn+1 )Ge(u n ) : e(u n )∂ ∆t ϕ n+1 . ( 83
)
Lemma 2. Let N ∈ N. Then we have

N -1 n=0 c( φn+1 )p n+1 -c( φn )p n p n+1 = 1 2 c( φN )p 2 N -c( φ0 )(p 0 ) 2 + 1 2 N -1 n=0 c( φn+1 ) p n+1 -p n 2 + 1 2 N -1 n=0 J ( φn+1 )(2p n+1 -p n )p n . ( 84
)
Lemma 3. Let N ∈ N. Then for the combination of the cross-terms, which are not multiplied by (α -1), we have

N -1 n=0 ( φn+1 ) 2 ∇p n+1 ∂ ∆t u n+1 + div ∂ ∆t (( φn+1 ) 2 u n+1 )p n+1 +D( φn+1 )∂ ∆t ϕ n+1 u n ∇p n+1 = N -1 n=0 div p n+1 ∂ ∆t (( φn+1 ) 2 u n+1 ) . ( 85
)
Lemma 4. Let N ∈ N. Then for the combination of the cross-terms, which are multiplied by (α -1), we have

N -1 n=0 ∇(( φn+1 ) 2 p n+1 )∂ ∆t u n+1 + p n+1 ∂ ∆t (( φn+1 ) 2 div u n+1 )- D( φn+1 )∂ ∆t ϕ n+1 p n+1 div u n = N -1 n=0 div p n+1 (( φn+1 ) 2 ∂ ∆t u n+1 ) . (86)
In the considerations which follow we will suppose that p bdry = 0. Working with the non-homogeneous pressure boundary condition would lead to long and cumbersome expressions but, on the other hand, their presence would not present any problem.

We introduce the incremental Lyapunov functional (or the incremental free energy functional) by

J N = 1 2 B A( φN )Ge(u N ) : e(u N ) + G c ε|∇ϕ N | 2 + 1 ε (1 -ϕ N ) 2 +c( φN )p 2 N dx - ∂ N B τ (t N ) • u N dS. ( 87
)
Then we have the following result Theorem 2. Let us suppose Hypothesis 1 and, in addition, p bdry = 0. Then the following estimate holds:

J N + ∆t N -1 n=0 B K ef f ( φn )∇p n+1 • ∇p n+1 dx ≤ J 0 + ∆t N -1 n=0 B K ef f ( φn ) G( φn ) • ∇p n+1 dx -∆t N -1 n=0 ∂ N B ∂ ∆t τ (t n+1 ) • u n+1 dS +∆t N -1 n=0 ∂ D B v inj (t n+1 )p n+1 dS + ∆t N -1 n=0 B q(t n+1 )p n+1 dx, ( 88
)
where J 0 is calculated using the initial values of the unknowns.

Proof. We consider system (78)-( 81) for times t n , n = 1, . . . , N. We test equation ( 78) by ∂ ∆t u n+1 , using equality (80) with ∂ ∆t ϕ n+1 and testing equation ( 81) by p n+1 . Then we sum up the obtained equalities from n = 0 to n = N -1 and, finally, sum up the three equalities obtained this way. Application of equalities ( 83)-( 86) yields the estimate (88).

Discretization and Solution Algorithm

For numerical treatment of the whole system, we work with Definition 1; namely equations ( 52)-( 54), in which the setting is time-discretized using a backward difference quotient. Then, we describe spatial discretization that is based on a Galerkin finite element scheme. Finally, we propose a combined Newton method that simultaneously solves for the nonlinear (forward) problem and the enforcement of the crack irreversibility.

Aspects of temporal discretization

All time derivatives are approximated with backward difference quotients as described in the theoretical part of the article. That is:

∂ t ϕ ≈ ϕ -Φ ∆t ,
where Φ := ϕ(t n-1 ) and the time step size is ∆t = t n -t n-1 . In a similar way the time derivative of the pressure equation is formulated. Moreover, we time-lag in terms of a linear extrapolation ϕ ex of the true ϕ in several terms in order to resolve numerical difficulties related to the indefiniteness of the Jacobian matrix. This is of most importance in the displacement equation. In [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF], it has been shown that this procedure is efficient and robust.

Spatial discretization

Spatial discretization is performed with conforming finite elements; namely the discrete spaces consists of bilinears using a mesh decomposed into quadrilaterals. This space is denoted as usually by Q c 1 . The discretization parameter is denoted by h and since the method is conforming, we have V h U ⊂ V U for the displacements, W h ⊂ W for the phase-field variable, and V h P ⊂ V P for the pressure. We employ a Galerkin scheme and formulate a semi-linear form

A(U h )(Ψ h ) with U h = (u h , ϕ h , p h ) ∈ V h U ×W h ×V h P
such that the discrete nonlinear system then reads: Definition 2 (of a single semi-linear form for the penalized formulation). In each time step t n , solve for a discrete solution

U h = (u h , ϕ h , p h ) ∈ V U × W × V P such that A(U h )(Ψ h ) = B A(ϕ h ex )Ge(u h ) : e(w h ) + (ϕ h ex ) 2 ∇p h + (α -1)∇((ϕ h ex ) 2 p h ) • w h dx - ∂ N B (τ • w h + p bdry w h n ) dS +Gc B - 1 ε (1 -ϕ h )ψ h + ε∇ϕ h • ∇ψ h dx + B δ(ϕ h -Φ)+ψ h dx + B D(ϕ h ex ) ((1 -α)p h div U + ∇p h • U ψ h dx + 1 2 B A(ϕ h ex )Ge(U) : e(U) -J (ϕ h ex )P (2p h -P ) ψ h + B c(ϕ h ex )p h + div ((ϕ h ex ) 2 u h ) + (α -1)(ϕ h ex ) 2 div u h π h dx +∆t B K ef f (Φ+)(∇p h -G(Φ+)) • ∇π h dx = ∆t B qπ h dx -∆t ∂ D B v inj πr dS + B c(Φ+)P + div (Φ 2 + U) + (α -1)Φ 2 + div U πr dx (89)
for all

Ψ h = (w h , ψ h , π h ) ∈ V h U × W h × V h P .
Remark 4 (Treatment of the crack irreversibility condition). The crack irreversibility can be either treated via penalization using the parameter δ as it has been used for related problems of pressurized and fluid-filled fractures [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF][START_REF] Wheeler | An augmented-Lagrangian method for the phase-field approach for pressurized fractures[END_REF]; or as recently proposed in [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF], employing a robust primal-dual active set strategy that can be interpreted as a semi-smooth Newton method [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF]. Further aspects are explained in more detail below in Section 6.4.

Newton's method for the nonlinear forward problem

In this section, we restrict our attention to the nonlinear solution of the main equations (i.e., the forward problem). We apply Newton's method to the fully-coupled problem in order to identify optimal convergence rates. Here, the Jacobian matrix is derived analytically in order to identify optimal Newton convergence whenever possible. The semi-linear form of the forward problem is defined in Definition 2 and neglects the penalization term. The nonlinear iteration reads:

Algorithm 1 (Newton iteration of the fully-coupled system). At a given time level; repeat the Newton iterations for k = 0, 1, 2, . . . :

1. Find dU h k := {du h k , dϕ h k , dp h k } by solving the linear system

A (U h k )(dU h k , Ψ h ) = -A(U h k )(Ψ h ) ∀Ψ ∈ V h × W h × V h . (90)
2. Find a step size 0 < ω ≤ 1 using line search (see Remark 5) to get

U h k+1 = U h k + ωdU h k , (91) with A
(U h k+1 )(Ψ h ) < A(U h k )(Ψ h
). Finish the Newton loop if the stopping criteria is fulfilled:

|A(U h k )(Ψ h )| < TOL .
In Algorithm 1, the Jacobian is derived by computing the Gâteaux derivative. The application to a semi-linear form reads:

A (U )(dU, Ψ) := lim s→0 1 s {A(U + sdU )(Ψ) -A(U )(Ψ)} = d ds A(U + sdU )(Ψ)| s=0 .
The Jacobian to Definition 2 reads:

A (U h )(dU h , Ψ h ) = B A(ϕ h ex )Ge(du h ) : e(w h ) + (ϕ h ex ) 2 ∇dp h + (α -1)∇((ϕ h ex ) 2 dp h ) • w h dx +G c B 1 ε (dϕ h )ψ h + ε∇dϕ h • ∇ψ h dx + B D(ϕ h ex ) ((1 -α)dp h div U + ∇dp h • U ψ h dx+ 1 2 B -J (ϕ h ex )P (2dp h ) ψ h + B c(ϕ h ex )dp h + div ((ϕ h ex ) 2 du h ) + (α -1)(ϕ h ex ) 2 div du h π h dx +∆t B K ef f (Φ + )(∇dp h -G(Φ + )) • ∇π h dx = 0 (92) for all Ψ h = (w h , ψ h , π h ) ∈ V h U × W h × V h P .
Remark 5 (Backtracking line search). A crucial role for (highly) nonlinear problems includes the appropriate determination of ω. A simple backtracking strategy is employed in the present paper and consists of modifying the update step in (91) as follows: For given ω ∈ (0, 1) (in our numerical tests, we choose ω = 0.6) determine the minimal l * ∈ N via l = 0, 1, . . . , N l , such that

R(U h k+1,l ) < R(U h k,l ), U h k+1,l = U h k+1 + ω l dU h k .
For the minimal l * , we set

U h k+1 := U h k+1,l * .
In this context, the nonlinear residual R(•) is defined as

R(U h ) := max i A(U h )(Ψ i ) -F (Ψ i ) ∀U h = (w h , ψ h , π h ) ∈ V h U × W h × V h P ,
where {Ψ i } denotes the nodal basis of V h U ×W h ×V h P . This algorithm works quite well for our problems and is applied to both the nonlinear forward model and the semi-smooth Newton method to realize the primal-dual active set.

A combined Newton method for the variational inequality system

We combine the Newton solver from the previous section with a semi-smooth Newton iteration for treating the crack irreversibility. Such an algorithm has been proposed in [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF]. Instead of penalizing the equations to solve the variational inequality, key idea of the primal-dual active set algorithm is that the domain is partitioned into an active set and an inactive set. In the active set, the constraint (here ϕ = Φ) is active, whereas in the inactive set we have ϕ < Φ and we solve the PDE. The interpretation as a semi-smooth Newton method allows detailed convergence analysis [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF] and fast convergence; i.e. a locally superlinear convergence.

With these preparations, we extend Algorithm 1 and define the combined Newton active-set iteration:

Algorithm 2 (Combined Newton of the fully-coupled system). At a given time level; repeat the Newton iterations for k = 0, 1, 2, . . . :

1. Assemble the residual A(U h k )(Ψ h ) with U h k = (u h k , ϕ h k , p h k )
formulated in Definition 2 consisting of displacements, phase-field, and pressure 2. Compute active set A k (see Remark 6 and [13]).

Assemble matrix A

(U h k )(•, Ψ h ) and right-hand side -A(U h k )(Ψ h )
4. Eliminate rows and columns in active set A k from matrix and right hand side to obtain

A (U h k )(•, Ψ h ) and -A(U h k )(Ψ h )
5. Solve the linear system using a GMRES method with block-diagonal preconditioning; namely find

dU h k such that A (U h k )(dU h k , Ψ h ) = -A(U h k )(Ψ h ) ∀Ψ ∈ V h × W h × V h . ( 93
)
6. Find a step size 0 < ω ≤ 1 using the line search algorithm explained in Remark 5 to compute

U h k+1 = U h k + ωdU h k , with A(U h k+1 )(Ψ h ) < A(U h k )(Ψ h
). Finish the Newton loop if both stopping criteria are fulfilled simultaneously:

A k+1 = A k and A(U h k )(Ψ h ) < TOL .
Remark 6 (Computing the active set). Determining the active set follows the procedure as it has been described, for example, in [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF]. On the discrete level, we check at each degree of freedom if the index i is in the active set A k :

(B -1 ) ii (R(U h k )) i + c(dU k h ) i > 0. ( 94
)
Otherwise, the index belongs to the inactive set I. Here, B is a mass matrix (details in [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF]), R(U h k ) the nonlinear residual from before, c > 0 a constant, and dU k h Newton's update. Each row and corresponding column of degree of freedoms i that are part of the active set are then eliminated from the Jacobian A (U h k )(•, Ψ h ) and the right hand side A(U h k )(Ψ h ) and thus the reduced systems are denoted by A (U h k )(dU h k , Ψ h ) and A(U h k )(Ψ h ), respectively. Inside Algorithm 2, we solve the linear equations with the generalized minimal residual (GMRES) method and block-diagonal preconditioning as implemented in [START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF] for the fully-coupled solution for {u, ϕ}. That program has been extended to solve the {u, ϕ, p} system. The (reduced) linear equation system (93) (AdU = b) has the following structure:

  A uu 0 B up C ϕu A ϕϕ M ϕp B pu 0 A pp     du h k dϕ h k dp h k   =   b u b ϕ b p   .
Here, the right hand side values {b u , b ϕ , b p } are assembled in terms of Newton's (reduced) residual A(U h k )(Ψ h ). The terms in the Jacobian matrix are specified through assembling the directional derivatives in (6.3). Let N u , N ϕ and N p be the dimensions of the spaces V h u , W h and V h p . Then,

A uu = (A(ϕ h ex )Ge(w j ) : e(w i )) Nu,Nu i,j=1 , (95)
B up = ((ϕ h ex ) 2 ∇π j + (α -1)∇((ϕ h ex ) 2 π j ) • w i ) Nu,Np i,j=1 , (96)
C ϕu = 0, (97)
A ϕϕ = G c ( 1 ε (ψ j )ψ i + ε∇ψ j • ∇ψ i ) Nϕ,Nϕ i,j=1 , (98)
M ϕp = (D(ϕ h ex )(((1 -α)π j div U + ∇π j • U)ψ i ) (99) -((J (ϕ h ex )P (2π j ))ψ i )) Nϕ,Np i,j=1 , (100)
B pu = (( div ((ϕ h ex ) 2 w j ) + (α -1)(ϕ h ex ) 2 div w j )π i ) Np,Nu i,j=1 , (101)
A pp = ((K ef f (Φ + )(∇π j -G(Φ + ))) • ∇π i ) Np,Np i,j=1 . (102)
The structure of the Jacobian A (U h k )(dU h k , Ψ h ) is important since with this information, a preconditioner P -1 can be constructed such that P -1 A has a moderate condition number. Here we simply choose

P -1 =   A -1 uu 0 0 0 A -1 ϕϕ 0 0 0 A -1 pp   .
This preconditioner is computed with an algebraic multigrid method from Trilinos [START_REF] Heroux | An overview of the trilinos project[END_REF] and yields a moderate number of iterations of the GMRES solver. For the examples in this paper we observed between 10 and 90 GMRES iterations per Newton step.

Numerical tests

In this final section, three different configurations of fluid-filled fractures in poroelasticity are considered. Specifically, Biot's coefficient is chosen as α = 1. All remaining parameters are summarized in Table 2. The domain is Ω = (0, 4) 2 m with boundary ∂Ω. As boundary conditions, we describe u = 0 on ∂Ω, ∂ n ϕ = 0 on ∂Ω and ∂ n p = 0 on ∂Ω (no flow). In order to study the robustness of our numerical solutions, we compute the first two examples on three different meshes using global mesh refinement. We emphasize that ε = 0.088m is fixed such that we are able to study h convergence.

SYMBOL
A point source injection is modeled as q = p ph d π exp(-d * (x -x 0 ) 2 ) with d = 10000 with x 0 = (2, 2). The fracture permeability is computed as K F ∼ w 2 /12η F where w is the width of the fracture; for a detailed derivation, we refer to [START_REF] Mikelić | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF]. The function c(ϕ + ) is a linear interpolation with c(0) = c F and c(1) = 1 M . As functionals, we observe the crack length and the maximum pressure over time.

A straight propagating fracture

In this first example, an initial crack is described with length l 0 = 0.4 on C = (1.8 -h, 2.2 + h) × (2h, 2 -h) ⊂ Ω. A single-phase fluid with source term q is injected into the middle of the domain in the point (2, 2).

Our findings can be observed in Figures 2345. The crack pattern is displayed in Figure 2. The pressure fields of the fracture and the reservoir interact. Specifically, the pressure at the tips of the fracture is very low and even becomes negative as observed in Figure 3. These findings are confirmed in Figure 4 in which the pressure is plotted over the middle cross section. This behavior has been already discussed in the literature by others [START_REF] Detournay | The near-tip region of a fluid-driven fracture propagating in a permeable elastic solid[END_REF][START_REF] Schrefler | On adaptive refinement techniques in multi-field problems including cohesive fracture[END_REF][START_REF] Secchi | A method for 3-D hydraulic fracturing simulation[END_REF]. The displacement field are illustrated in Figure 5. We observe the typical shape of the crack opening displacement; namely zero opening at the tips of the fracture and largest opening in the middle. This phenomenon was predicted by the theoretical calculations made in [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF] in their Section 2.4 for 2D configurations. Finally, we observe in Figure 6 the crack length and the total pressure on three different meshes in order to study the robustness of our solutions. Our findings indicate that our method is computationally stable for h-refinement while keeping ε fixed. The units of horizontal and vertical axes are m and P a, respectively. Here, we observe large pressures in the fracture and low (negative) pressures at the (propagating) tips. As long as the fracture is not propagating, i.e., between T = 1s and T = 10s, the pressure increases. While the fracture is propagating we observe a slight pressure decrease. In addition, while the fracture is propagating negative pressures at the tips remain but they do not decrease further. For fixed ε, we observe h-convergence in the length as well as in the pressure. In particular, we observe a pressure drop while the crack is growing. On the coarse mesh the pressure again increases while the pressure is monotonically decreasing on the two fine meshes.

Propagation of multiple fractures in a homogeneous porous medium

In this second test, the domain, boundary conditions and the first fracture (the horizontal one) remain the same as in the first example. We now add a second, vertically-aligned fracture on (2.6 ± h, 2.0 ± 0.2)m. We only do inject a fluid into fracture 1 (like in the first example). Consequently, the second fracture can be seen as a natural non-pressurized fracture. Analyzing our findings, the following conclusions can be inferred: In Figure 7, we observe the crack pattern evolution and see that only fracture 1 starts propagating because fracture 2 is not pressurized. They join at time T = 7. Next, the horizontal crack continues further propagating. Later, the second crack also starts propagating. That fracture 2 is not pressurized can be justified in Figure 8 since high pressure is in the first time steps only being seen in fracture 1. When they join, also fracture 2 is pressurized and fluid flows inside. Similarly to Example 7.1, the pressure field is negative at the tips of the fractures. Finally, convergence studies of the half crack length of the horizontal fracture and maximum pressure versus time on three different mesh levels are shown in Figure 9. For fixed ε, we observe h-convergence in the length as well as in the pressure. In particular, we observe a huge pressure drop at T = 7s when the horizontal crack joins the vertical crack. In this example, we also observe that a too coarse mesh does might give wrong answers as here seen for the pressure. The pressure increase towards T = 50 is due to the fact that this crack reaches the boundary and lead to disputable results.

Fracture propagation in a heterogeneous porous medium

In this final example, the reservoir permeability K R varies. This is an important example motivated by practical applications. Examples with varying geomechanics parameters have been presented in [START_REF] Mikelić | A quasistatic phase field approach to pressurized fractures[END_REF][START_REF] Heister | A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach[END_REF] as well as in [START_REF] Wick | Pressurized-fracture propagation using a phase-field approach coupled to a reservoir simulator[END_REF] to couple phase-field fracture propagation with a reservoir simulator.

Except of the reservoir permeability, the configuration, boundary conditions, model and material parameters are the same as in Example 7.1. Analyzing our results, we observe that the fracture now grows on a non-planar way due to the non-constant permeability field (see Figure 11). The pressure field is displayed in Figure 12. Again, low pressures occur at the tips of the fracture. We note that the displacements are zero at the tips of the fracture and the highest value is obtained in the middle. As in the first example, this corresponds to the theoretical calculations of [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF] in their Section 2.4.

Conclusion

In this paper, we formulated a fully-coupled phase-field model describing propagating cracks in a poroelastic medium. The model consists of coupled nonlinear PDEs for the displacement, the pressure and the phase field unknowns. The phase-field algorithm is based on an incremental formulation and existence of a weak solution is established. Furthermore, the corresponding Lyapunov functional, based on the free energy, is constructed. Numerical tests demonstrate the correctness of the model. Specifically, this approach can treat crack growth in heterogeneous porous media, pressurized crack evolutions and multiple interacting cracks allowing them to join and branch. Moreover, the numerical simulations confirm findings from literature that low (negative) pressures arise at fracture tips.

  is denoted by ∂B, which is divided into open 2D surfaces ∂ D B and ∂ N B, with smooth boundaries, in the way that ∂B = ∂ D B ∪ ∂ N B. On ∂ D B Dirichlet boundary conditions will be imposed and on ∂ N B Neumann conditions, respectively. We assume that meas(∂ D B) > 0. Boundary conditions on ∂ D B for the Biot system involve displacements and the normal pressure fluxes and on ∂ N Ω traction as well as pressure.

Proposition 1 .

 1 Let us suppose smooth data and a smooth solution {u, ϕ, p} for problem (32)-(38). Then the following energy equality holds + )Ge(u) : e(u)

Definition 1 (

 1 of a penalized approximation). Let us suppose that the data satisfy Hypothesis 1, a penalization parameter δ ∈ R and in particular, let δ := M ∈ N in this section. Let φ = inf{1, ϕ + } and let V M U = span {w r } r=1,...,M , V M P = span {π r } r=1,...,M and W M = span {ψ r } r=1,...,M . The triple {u M , ϕ M , p M }, u M = M r=1 a r w r , p M = p bdry + M r=1 d r π r and ϕ M = Φ + M r=1 b r ψ r , is a finite dimensional approximative solution for problem (42)-(50) if it satisfies the discrete variational Proof. (of Proposition 2) Let ξ = {a r , b r , d r } r=1,.

Figure 2 :

 2 Figure 2: Example 1: Crack pattern at times T = 1, 10, 20, 30s using the phase-field variable ϕ. In the crack (red colored) region ϕ = 0 and in the unbroken material, ϕ = 1. In between 0 < ϕ < 1, which denotes the smooth transition zone.

Figure 3 :

 3 Figure 3: Example 1: Pressure (in P a) at times T = 1, 10, 30, 40s. Here, we observe a large pressure difference with negative pressures in front of the tips of the fracture.

Figure 4 :

 4 Figure 4: Example 1: Pressure plot over the cross section {(x, 2)|0 ≤ x ≤ 4} at times T = 1, 10, 30, 40s.The units of horizontal and vertical axes are m and P a, respectively. Here, we observe large pressures in the fracture and low (negative) pressures at the (propagating) tips. As long as the fracture is not propagating, i.e., between T = 1s and T = 10s, the pressure increases. While the fracture is propagating we observe a slight pressure decrease. In addition, while the fracture is propagating negative pressures at the tips remain but they do not decrease further.

Figure 5 :

 5 Figure 5: Example 1: y-displacement field (in m) at times T = 1, 10, 30, 40s. Here, the displacements are zero at the tips of the fracture and the highest value is obtained in the middle. As in the first example, this corresponds to the theoretical calculations of [20] in their Section 2.4.

Figure 6 :

 6 Figure 6: Example 1: Convergence studies of crack length and maximum pressure versus time on three different mesh levels.For fixed ε, we observe h-convergence in the length as well as in the pressure. In particular, we observe a pressure drop while the crack is growing. On the coarse mesh the pressure again increases while the pressure is monotonically decreasing on the two fine meshes.

Figure 7 :

 7 Figure 7: Example 2: Crack pattern at times T = 1, 6, 7, 20, 30, 50s with two fractures. Fracture 1 (the horizontal fracture) joins fracture 2 (the vertical fracture) at T = 7s. The crack pattern is displayed with the help of the phase-field variable ϕ. In the crack (red colored) region ϕ = 0 and in the unbroken material, ϕ = 1. In between 0 < ϕ < 1, which denotes the smooth transition zone.

Figure 8 :

 8 Figure 8: Example 2: Pressure (in P a) at times T = 1, 6, 7, 20, 30, 50s. Again, we observe a large pressure difference with low pressures (in blue color) in front of the fracture tips. In fact, the pressure is lowest first at right when fracture 1 joins fracture 2. Then, fracture 1 further grows and we observe lowest pressure at T = 20s at left in front of the growing fracture.Later, the two fracture tips of the second fracture start growing and now the lowest pressures are seen here at T = 50s. In this final figure properly the pressure might be still very low for the fracture 1 tip at left, but this fracture is leaving the domain in this moment. Here, we notice that boundary interactions/influences are still subject of current research. But in the last figure, our primary focus is on the tips of fracture 2.

Figure 9 :

 9 Figure 9: Example 2: Convergence studies of half crack length of the horizontal fracture and maximumpressure versus time three different mesh levels. For fixed ε, we observe h-convergence in the length as well as in the pressure. In particular, we observe a huge pressure drop at T = 7s when the horizontal crack joins the vertical crack. In this example, we also observe that a too coarse mesh does might give wrong answers as here seen for the pressure. The pressure increase towards T = 50 is due to the fact that this crack reaches the boundary and lead to disputable results.

Figure 10 :

 10 Figure 10: Example 3: Initial crack, pressure and permeability distribution

Figure 11 :

 11 Figure 11: Example 3: Crack pattern at times T = 10, 20, 30, 50s using the phase-field variable ϕ. In the crack (red colored) region ϕ = 0 and in the unbroken material, ϕ = 1. In between 0 < ϕ < 1, which denotes the smooth transition zone.

Figure 13 :

 13 Figure 13: Example 3: y-displacement field (in m) at times T = 10, 20, 30, 50s.

Table 1 :

 1 Unknowns and effective coefficients

	SYMBOL	QUANTITY	UNITY
	u	displacement	m
	p	fluid pressure	Pa
	σ por	total poroelasticity tensor	Pa
	e(u) = (∇u + ∇ τ u)/2 linearized strain tensor	dimensionless
	KR	permeability of the poroelastic domain Darcy
	KF	fracture permeability	Darcy
	α	Biot's coefficient	dimensionless
	ρ b	bulk density	kg/m 3
	ηF	fracture fluid viscosity	kg/m sec
	ηR	pore fluid viscosity	kg/m sec
	M	Biot's modulus	Pa
	cF	fracture fluid compressibility	1/Pa
	G	Gassman rank-4 tensor	Pa
	ρF	reference state fracture fluid density	kg/m 3
	ρR	reference state pore fluid density	kg/m 3
	g	gravity	m/sec 2

•) the trace operator, and µ and λ are Lamé's parameters, which are linked to Young's modulus E and Poisson's ratio ν s by µ = E/(2(1 + ν s )) and λ = Eν s / (1 + ν s )(1 + 2ν s ) .

Table 2 :

 2 Parameters for the numerical examples

Lemma 1.4. (see e.g. R. Temam, Navier-Stokes Equations, page 164) Let X be a finite dimensional Hilbert space with scalar product (•, •)X and norm || • ||X and let P be a continuous from X into itself such that (P(ξ), ξ)X > 0 for ||ξ||X = R > 0.Then there exists ξ ∈ X , ||ξ||X ≤ R, such that P(ξ) = 0.

* The research by A. Mikelic, M. F. Wheeler and T. Wick was partially supported by ConocoPhillips grant UTA13-001170 AMD 1 and Statoil grant UTA13-000884: WR DTD 2.13.14. A. Mikelic and T. Wick would like to thank Institute for Computational Engineering and Science (ICES), UT Austin for hospitality during their visits in February, June and August 2015 and support through a JT Oden fellowship. M. F. Wheeler was also supported by Aramco grant UTA11-000320; 1ST and T. Wick by the Austrian Academy of Sciences.

with

Proposition 2. We suppose Hypothesis 1. Then there exists a finite dimensional penalized approximative solution for problem (52)-(54) that satisfies the a priori estimate

where C is independent of M and δ.

After inserting (61)-( 62) into (60), we get

Next we have

Estimating different the terms is straightforward:

Therefore, we have

It remains to consider the third component P 3 of P :

The idea is to put (63), (66) and (67) together and to obtain a convenient lower bound. A quick inspection of the obtained estimates shows that it remains to control the cross terms in (P 1 (ξ), ξ 1 )


Metod Saniga
email: msaniga@astro.sk

Combinatorial Intricacies of Labeled Fano Planes

Keywords:

are published or not. The documents may come   L'archive ouverte pluridisciplinaire

It is well known [START_REF] Hall | On identifying PG(3,2) and the complete 3-design on seven points[END_REF][START_REF] Lindner | Steiner triple systems having a prescribed number of triples in common[END_REF][START_REF] Polster | Yea why try her raw wet hat: A tour of the smallest projective space[END_REF] that there are thirty different Fano planes on a given seven-element set. Slightly rephrased, there are thirty different ways to label the points of the Fano plane by integers from 1 to 7, two labeled Fano planes having zero, one or three lines in common and each line occurring in six Fano planes. The set of thirty labeled Fano planes can be uniquely partitioned into two sets, say A and B, of fifteen elements each, such that any two labeled Fanos in the same set have just one line in common. Following the labeling adopted by Polster [START_REF] Polster | A Geometrical Picture Book[END_REF], the two sets read where for triples of integers we use a shorthand notation {a, b, c} = abc. Given a line abc, where, without loss of generality, we can take 1 ≤ a < b < c ≤ 7, we shall distinguish between the cases when a + b = c and a + b = c and call the former/latter ordinary/defective [START_REF] Saniga | Cayley-Dickson algebras and finite geometry[END_REF]; for the sake of convenience, in the above-given sets A and B all ordinary lines are denoted by overbars. Furthermore, a point of the labeled Fano plane is said to be of order s, 0 ≤ s ≤ 3, if there are s defective lines passing through it; hence, in addition to two different kinds of lines, a labeled Fano plane can potentially feature up to four distinct types of points. A detailed inspection of each of the 30 labeled Fano planes given above shows that they fall, in terms of such structural refinement, into eight different types as portrayed in Figure 1 and summarized in Table 1. It is also interesting to compare the distributions of types within each 15-element set, given in Tables 2 and3, as well as cardinalities of individual types, listed in Table 4; note a pronounced asymmetry between sets A and B.

A 1 : {124,
Figure 1: A diagrammatic illustration of representatives of eight distinct types of numbered Fano planes. A point of order three, two, one and/or zero is represented by red, blue, green and/or yellow color, respectively; heavy lines are defective. The types α and γ do not exist, i. e. it is impossible label the Fano plane in such a way that all or just five of its lines would be ordinary.

Table 1: The eight distinct types of labeled Fano planes, each being uniquely characterized by the number of points of every particular order. To understand better the symmetry pattern, we arrange complementary types into pairs and list also the two non-existent types α and γ.

Points of Order

Type 0 1 2 3 (α) (7) (0) (0) (0) α ′ 0 0 0 7 β 4 3 0 0 β ′ 0 0 3 4 (γ) (2) (4) (1) (0) γ ′ 0 1 4 2 δ 1 1 3 3 0 δ ′ 1 0 3 3 1 δ 2 0 6 0 1 δ ′ 2 1 0 6 0
γ ′ β ′ β ′ γ ′ β ′ γ ′ γ ′ β ′ γ ′ β β ′ α ′ β ′ δ 2 β ′
α ′ α ′ γ ′ δ ′ 2 α ′ δ ′ 2 α ′ γ ′ γ ′ δ ′ 1 γ ′ δ 1 γ ′ α ′ Table 4: Cardinalities of individual types of Fano planes. Type (α) β (γ) δ 1 δ 2 α ′ β ′ γ ′ δ ′ 1 δ ′ 2 Set A - 1 - 0 1 1 7 5 0 0 Set B - 0 - 2 0 5 0 5 1 2 Total - 1 - 2 1 6 7 10 1 2
In addition to classifying planes, we can also classify lines in terms of the types of six planes passing through each of them, the corresponding findings being given in Table 5. We find that nine ordinary lines are of five different kinds, whereas the remaining 26 defective lines yield as many as ten distinct types. We further observe that no ordinary line exhibits a plane of type α ′ and there is only one defective line that features plane of type β. It is also worth pointing out that it is only one plane in set A that is devoid of ordinary lines, whereas in set B one finds four such planes.

To conclude the paper, let us have a look at Figure 1 and focus on that sub-geometry in each labeled Fano plane representative which is formed by defective lines and points lying on them. We shall find the Fano plane itself and, using the language of [START_REF] Grannell | Three-line chromatic indices of Steiner triple systems[END_REF], all its derivatives depicted in Figure 2 that play a crucial role in classifying Steiner triple systems (see, for example, [START_REF] Horak | Counting frequencies of configurations in Steiner triple systems[END_REF][START_REF] Forbes | Configurations and trades in Steiner triple systems[END_REF]); as readily discerned from Figure 1, the semihead lies in a type-β ′ plane, the mia in a type-γ ′ plane, a sail is hosted by the type-δ ′ 1 plane, the Pasch configuration is located in a type-δ ′ 2 plane, the triangle in a type-δ 1 plane and the Fano plane itself coincides with a type-α ′ plane. It should also be pointed out that if one instead considers sub-geometries formed by ordinary lines, then both the Fano plane and its derivate mia are absent.

Obviously, any finite geometry (point-line incidence structure) with lines of size three can be looked at this way. We have already performed the corresponding examination of the Möbius-Kantor (8 3 )-configuration, the Pappus (9 3 )-configuration as well as the Desargues (10 3 )-configuration and plan to publish its outcome in a separate paper.

  136, 157, 235, 267, 347, 456}, B

	A 6 : {124, 137, 156, 236, 257, 345, 467},	B 6 : {125, 137, 146, 236, 247, 345, 567},
	A 7 : {123, 147, 156, 245, 267, 346, 357},	B 7 : {123, 145, 167, 247, 256, 346, 357},
	A 8 : {124, 135, 167, 237, 256, 346, 457},	B 8 : {126, 135, 147, 237, 245, 346, 567},
	A 9 : {126, 137, 145, 235, 247, 346, 567},	B 9 : {124, 137, 156, 235, 267, 346, 457},
	A 10 : {123, 145, 167, 246, 257, 347, 356},	B 10 : {123, 146, 157, 245, 267, 347, 356},
	A 11 : {126, 134, 157, 237, 245, 356, 467},	B 11 : {125, 134, 167, 237, 246, 356, 457},
	A 12 : {125, 137, 146, 234, 267, 356, 457},	B 12 : {126, 137, 145, 234, 257, 356, 467},
	A 13 : {123, 146, 157, 247, 256, 345, 367},	B 13 : {123, 147, 156, 246, 257, 345, 367},
	A 14 : {127, 134, 156, 235, 246, 367, 457},	B 14 : {126, 134, 157, 235, 247, 367, 456},
	A 15 : {126, 135, 147, 234, 257, 367, 456},	B 15 : {127, 135, 146, 234, 256, 367, 457},

1 : {127, 136, 145, 235, 246, 347, 567}, A 2 : {127, 136, 145, 234, 256, 357, 467}, B 2 : {125, 136, 147, 234, 267, 357, 456}, A 3 : {125, 136, 147, 237, 246, 345, 567}, B 3 : {124, 136, 157, 237, 256, 345, 467}, A 4 : {125, 134, 167, 236, 247, 357, 456}, B 4 : {127, 134, 156, 236, 245, 357, 467}, A 5 : {127, 135, 146, 236, 245, 347, 567}, B 5 : {124, 135, 167, 236, 257, 347, 456},

Table 2 :

 2 Types of Fano planes in set A.

	Plane 1	2	3	4	5	6	7	8	9 10 11 12 13 14 15
	Type

Table 3 :

 3 Types of Fano planes in set B.

	Plane 1	2	3	4	5	6	7	8	9 10 11 12 13 14 15
	Type δ 1

Acknowledgment

This work was partially supported by the VEGA Grant Agency, Project 2/0003/13, as well as by the Austrian Science Fund (Fonds zur Förderung der Wissenschaftlichen Forschung (FWF)), Research Project M1564-N27. I am grateful to my friend Petr Pracna for electronic versions of the figures.

Pasch Mia


A Meziane
email: meziane.ahmed@univ-oran.dz

M Tahiri

Geometry and off-shell nilpotency for N = 1 supersymmetric

Keywords: Global supersymmetry, Auxiliary fields, Off-shell nilpotent BRST and anti-BRST algebra, Quantum action PACS numbers: 11.15.-q, 11.30.Pb, 12.60.Jv

We show that for N = 1 supersymmetric Yang-Mills theory it is possible to build an off-shell nilpotent BRST and anti-BRST algebra in terms of a BRST superspace formalism. This is based on the introduction of the basic fields of the quantized theory together with an auxiliary real field via the lowest components of the superfield components of a superYang-Mills connection. Here, the associated supercurvature is constrained by horizontality conditions as in ordinary Yang-Mills theory. We also show how the off-shell BRST-invariant quantum action can be constructed starting from a gauge-fixed superaction.

I. INTRODUCTION

It is well known that the quantization of globally supersymmetric gauge theories has been considerably studied a long time ago. Several approaches have been proposed with various methods in order to perform the quantization of such theories (for a review see e.g. Ref. [START_REF] Henneaux | Quantization of Gauge Systems[END_REF] and references therein).

In the component field formalism, the supersymmetry algebra achieved without auxiliary fields closes only on-shell. This can be explained by the fact that the supersymmetry transformations of the models are realized nonlinearly and therefore the main problem affecting such theories is linked to its algebraic structure which involves equations of motion and field dependent gauge transformations. This gives rise to an infinite dimensional algebra, even if auxiliary fields can be introduced to put the formalism off-shell [START_REF] Breitenlohner | and N = 2 Su-persymmetric Yang-Mills theories in the Wess-Zumino gauge[END_REF]. To avoid these difficulties it has been proposed in Ref. [START_REF] Maggiore | [END_REF] the construction of a generalized BRST operator by collecting together all the symmetries forming the theory, namely ordinary BRST, supersymmetries and translations. According to this procedure the role of the auxiliary fields is covered by the external sources coupled to the nonlinear variations of the quantum fields. This approach has been already successfully applied to supersymmetric [4,5], ordinary [6] gauge field theories as well as to non gauge field theories [7]. Let us mention that in the Wess-Zumino model it is only with auxiliary fields that one can obtain a tensor calculus (for a review see e.g. Ref. [8] and references therein). . Another possibility to solve the problem of the quantization of theories with on-shell algebra is to introduce the Batalin-Vilkovisky (BV) formalism [9]. The BV formalism is a very general covariant Lagrangian approach which overcomes the need of closed classical algebra by a suitable construction of BRST operator. The construction is realized by introducing a set of the so-called antifields besides the fields occurring in the theory. The elimination of these antifields at the quantum level via a gauge-fixing procedure leads to the quantum theory in which effective BRST transformations are nilpotent only on-shell. Let us note, that the BV approach can be used to obtain the on-shell BRST invariant gauge fixed action for N = 1 supersymmetric Yang-Mills theory in four dimensions without requiring the set of auxiliary fields [10].

Another interesting approach with infinite number of auxiliary fields has been developed in the context of harmonic superspace [11]. In this framework quantization of supersym-metric theories has been discussed for various supersymmetries [START_REF] Galperin | Harmonic Superspace[END_REF].

On the other hand, it is also known that the extension of spacetime with two ordinary anticommuting coordinates to a (4, 2)-dimensional superspace [START_REF] Bonora | [END_REF] leads in Yang-Mills type theories to incorporate the gauge fields, the ghost and anti-ghost fields into a natural gauge superconnection, see also Ref. [14] and references therein. In such a superspace formalism the BRST and anti-BRST transformations are derived systematically from the horizontality conditions imposed on the supercurvature.

Let us note that the natural geometrical way to derive the BRST structure of general gauge theories is to work, in the same spirit as in Yang-Mills type theories, by using the superconnection formalism. Within this framework and in contrast to what is done in Yang-Mills theories, all superfield components of the supercurvature cannot be constrained through horizontality conditions. This is a consequence of the fact that the gauge theories we consider are reducible and/ or open. It is the consistency of the Bianchi identities which is guaranteed by the remaining superfield components of the supercurvature. Their lowest components allow the introduction of auxiliary fields. These together with the fields given by the lowest components of the superfield components of the superconnection represent the basic fields of the quantized gauge theory. The off-shell nilpotency of the BRST and anti-BRST transformations of these fields is automatically ensured, thanks to the structure equations and the Bianchi identities. The BRST and anti-BRST operators are related as usual to the partial derivatives with respect to the two anti-commuting coordinates of the superspace. Essentially, the introduction of the auxiliary fields gives rise to the construction of the off-shell BRST invariant quantum action. As shown in Ref. [15] for the case of non-Abelian BF theory where the classical gauge algebra is reducible and Ref. [16] for the case of the simple supergravity where the classical gauge algebra is open, the superspace formalism has been used in order to realize the BRST structure of such theories. It leads to recast all the fields in geometrical way and to introduce auxiliary fields ensuring the off-shell invariance of the quantum action.

Our main aim in this paper consists to apply this formalism for discussing the off-shell nilpotent version of the BRST and anti-BRST transformations for global N = 1, 4D supersymmetric Yang-Mills theory where the classical gauge algebra is open [8,[START_REF] Wess | Supersymmetry and Supergravity[END_REF]. Let us mention that in Ref. [15,16] the superspace formalism has been applied successfully to theories with local symmetry while in the present work we are interested to apply this formalism to a global supersymmetry. The classical action for the N = 1 supersymmetric Yang-Mills in four dimensions is given by [START_REF] Wess | Supersymmetry and Supergravity[END_REF] S 0 = dx 4 tr( -

1 4 F µν F µν -λγ µ D µ λ), (1)
where tr denotes the trace over the gauge algebra, (A µ , λ) is the gauge multiplet, the field ] is the covariant derivative. By construction, in addition to the ordinary Yang-Mills symmetry, the action ( 1) is invariant under the supersymmetry transformations

strenght is F µν = ∂ µ A ν -∂ ν A µ + [A µ , A ν ] and D µ = ∂ µ + [A µ ,
δ ξ A µ = iξγ µ λ + iξγ µ λ, δ ξ λ = σ µν F µν ξ, (2)
where σ µν ≡ 1 2 [γ µ , γ ν ] and ξ is a spin 1/2 valued infinitesimal supersymmetry parameter. In the following we shall call the superspace obtained by enlarging spacetime with two ordinary anticommuting coordinates BRST superspace, in order to distinguish it from superspace of supersymmetric theories. Let us recall that a full off-shell structure of any supersymmetric field theory most naturally exhibits itself in superspace, provided the superfield formulation of the theory in terms of unconstrained superfields is available. The BRST superspace formalism we present here permits us to derive the off-shell nilpotent BRST and anti-BRST algebra of quantized N = 1, 4D supersymmetric Yang-Mills theory.

In particular, it gives another possibility leading to the minimal set of auxiliary fields in such theory.

Our paper is organized as follows: In Section II, the BRST superspace approach and horizontality conditions for N = 1, 4D supersymmetric Yang-Mills theory are discussed.

We also show how the various fields of such theory and their off-shell nilpotent BRST and anti-BRST transformations can be determined via a BRST superspace formalism. The construction of the BRST-invariant quantum action for N = 1 super Yang-Mills theory in terms of this off-shell structure is described in Section III. Section IV is devoted to concluding remarks.

II. OFF-SHELL NILPOTENT BRST ALGEBRA

Let Φ be a super Yang-Mills connection on a (4, 2)-dimensional BRST superspace with coordinates z=(z M )=(x µ , θ α ), where (x µ ) µ =1,..,4 are the coordinates of the spacetime manifold and (θ α ) α =1,2 are ordinary anticommuting variables. Acting the exterior covariant superdifferential D on Φ we obtain the supercurvature Ω satisfying the structure equation, Ω = dΦ+(1/2)[Φ,Φ], and the Bianchi identity, dΩ+[Φ,Ω]= 0. The superconnection Φ as 1-superform on the BRST superspace can be written as

Φ = dz M (Φ i M I i + Φ µ M P µ + Φ a M Q a ), (3)
where {I i } i=1,..,d=dimG and {P µ , Q a } µ=1,...4;a=1,..,4 are the generators of the internal symmetry group (G) and the N = 1 supersymmetric group (SG) respectively. They satisfy the following commutation relations

[I i , I j ] = f k ij I k , [I i , P µ ] = [Q a , P µ ] = [P µ , P ν ] = 0, [Q a , Q b ] = 2(γ µ ) ab P µ , [I i , Q a ] = b * i Q a , (4)
where {γ µ } µ=1,.. [START_REF] Wess | Supersymmetry and Supergravity[END_REF][START_REF] Weinberg | The Quantum Theory of Fields[END_REF]. Note that the Grassmann degrees of the superfield components of Φ are given by

|Φ i M |=|Φ µ M | = m, |Φ a M | = m + 1 (mod 2), where m = |z M | (m = 0 for M = µ and m = 1 for M = α), since Φ is an even 1-superform.
However, we assign to the anticommuting coordinates θ 1 and θ 2 the ghost numbers (-1) and (+1) respectively, and ghost number zero for an even quantity: either a coordinate, a superform or a generator. These rules permit us to determine the ghost numbers of the superfields

( Φ i µ , Φ ν µ , Φ a µ , Φ i 1 , Φ i 2 , Φ ν 1 , Φ ν 2 , Φ a 1 , Φ a
2 ) which are given by (0, 0, 0, +1, -1, +1, -1, +1, -1)

Upon expressing the supercurvature Ω as

Ω = 1 2 dz N ∧ dz M Ω M N = 1 2 dz N ∧ dz M (Ω i M N I i + Ω µ M N P µ + Ω a M N Q a ), (5)
we find from the structure equation

Ω µν = ∂ µ Φ ν -∂ ν Φ µ + [Φ µ , Φ ν ], (6a)
Ω µα = ∂ µ Φ α -∂ α Φ µ + [Φ µ , Φ α ], (6b)
Ω αβ = ∂ α Φ β + ∂ β Φ α + [Φ α , Φ β ]. (6c)
Similarly, the Bianchi identity becomes

D µ Ω νκ + D κ Ω µν + D ν Ω κµ = 0, (7a)
D α Ω µν -D ν Ω µα + D µ Ω να = 0, (7b)
D α Ω βγ + D β Ω αγ + D γ Ω αβ = 0, (7c)
D µ Ω αβ -D α Ω µβ -D β Ω µα = 0, (7d)
where

D M = ∂ M + [Φ M , .
] is the M covariant superderivative. Now, we shall search for the constraints to the supercurvature Ω in which the consistency with the Bianchi identities (7) is ensured. This requirement ensures then the off-shell nilpotency of the BRST and anti-BRST algebra. The full set of supercurvature constraints turns out to be given by

Ω µα = 0, Ω αβ = 0. ( 8
)
It is easy to check the consistency of this set of supercurvature constraints through an analysis of the Bianchi identities. Indeed, we remark that identities (7c) and (7d) are automatically satisfied because of the constraints (8) while identities (7a) and (7b) yield further restrictions on supercurvature Ω

Ω κ µν = 0, Ω a µν = 0. (9)
At this point, let us mention that the consistency of the horizontability conditions ( 8) and ( 9) with the Bianchi identities ( 7), as we will see later, guarantees automatically the off-shell nilpotency of the BRST and anti-BRST transformations on all the fields belonging to N = 1 super Yang-Mills theory.

Now in order to derive the off-shell BRST structure of N = 1 super Yang-Mills theory using the above BRST superspace formalism, it is necessary to give the geometrical interpretation of the fields occurring in the quantization of such theory. Besides the gauge potential Φ i µ| = A i µ , there exists the superpartener Φ i a| = λ i a of A µ and an auxiliary real field Φ i | = Λ i which are introduced via the field redefinitions

Φ ai = -1 4 (γ µ ) ab [Q b , Φ i µ ] and Φ i = 1 4 δ a b [Q a ,
Φ bi ] respectively. The components Φ ν µ| and Φ a µ| are identified to zero. Note that the symbol "|" indicates that the superfield is evaluated at θ α = 0. Furthermore, we introduce the following now:

Φ i 1| = c i 1 is the ghost for Yang-Mills sym- metry, Φ i 2| = c i 2 is the antighost of c i 1 , B i = ∂ 1 Φ i 2| is the associated auxiliary field, Φ a 1| = χ a 1 is the supersymmetric ghost, Φ a 2| = χ a 2 , is the antighost of χ a 1 , G a = ∂ 1 Φ a 2| is the associated auxiliary field, Φ µ 1| = ξ µ 1 is the translation symmetry ghost, Φ µ 2| = ξ µ 2 , is the antighost of ξ µ 1 and E µ = ∂ 1 Φ µ 2|
is the associated auxiliary field. Let us mention that the symmetry ghosts and antighosts χ a α are commuting fields while the others c i α and ξ µ α are anticommuting. The action of the N = 1 supersymmetric generators {P µ , Q a } on these fields is given by

[P µ , X] = ∂ µ X, [Q a , A i µ ] = -(γ µ ) ab λ bi , [Q a , λ bi ] = - 1 2 (σ µν ) b a F i µν + δ b a Λ i , [Q a , c i α ] = (γ µ ) ab χ b α A i µ , [Q a , Λ i ] = (γ µ ) ab D µ λ bi , [Q a , F i µν ] = (γ µ ) ab (D ν λ) bi -(γ ν ) ab (D µ λ) bi , ( 10
)
where X is any field of the theory.

It is worthwhile to mention that we are interested in our present investigation on the global supersymmetric transformations, so that the parameters of the N = 1 supersymmetric and translation groups must be space-time constant, i.e.

∂ µ χ a α = 0, ∂ µ ξ ν α = 0. ( 11
)
Using the above identifications with (11) and inserting the constraints ( 8) and ( 9) into the Eqs. ( 6) and (7b), we obtain

∂ α Φ i µ| = (D µ c α ) i -ξ ν α [P ν , A i µ ] -χ a α [Q a , A i µ ], ∂ α Φ i β| + ∂ β Φ i α| = -[c α , c β ] i -ξ ν α [P ν , c i β ] -χ a α [Q a , c i β ] -ξ ν β [P ν , c i α ] -χ a β [Q a , c i α ], ∂ α Φ a β| + ∂ β Φ a α| = 0, ∂ α Φ µ β| + ∂ β Φ µ α| = -2χ a α (γ µ ) ab χ b β , ∂ α Ω i µν| = -[c α , F i µν ] -ξ τ α [P τ , F i µν ] -χ a α [Q a , F i µν ], ∂ α Φ ai | = 1 4 (γ µ ) ab [Q b , ∂ α Φ i µ| ], ∂ α Φ i | = - 1 4 (δ a b )[Q a , ∂ α Φ bi | ]. ( 12
)
We also realize the usual identifications:

Q α (X | ) = ∂ α X |
, where X is any superfield and

Q = Q 1 (Q = Q 2 )
is the BRST (anti-BRST) operator. Inserting Eq. ( 10) into ( 12), and evaluating these at θ α = 0, we find the following BRST transformations

QA i µ = (D µ c) i -ξ ν ∂ ν A i µ + χγ µ λ i , Qλ i a = -f i jk c j λ k a -ξ µ ∂ µ λ i a + χ a Λ i + 1 2 (χσ µν ) a F i µν , Qc i = - 1 2 f i jk c j c k -ξ µ ∂ µ c i + χγ u χA i µ , Qξ µ = -χγ u χ, QΛ i = -f i jk c k Λ j -ξ ρ ∂ ρ Λ i -χγ µ D µ λ i , QF i µν = -f i jk c k F j µν -ξ ρ ∂ ρ F i µν -χ a {(γ µ ) ab (D ν λ) bi -(γ ν ) ab (D µ λ) bi }, Qχ a = 0, Qc i = B i , QB i = 0, Qξ µ = E µ , QE µ = 0, Qχ a = G a , QG a = 0, ( 13
)
and also the anti-BRST transformations, which can be derived from ( 13) by the following mirror symmetry of the ghost numbers given by :

X → X if X = A i µ , λ i a , Λ i ; X → X if X = Q, c i , B i , ξ µ , E µ , χ a , G a and X = X where B i + B i = -f i jk c k c j -ξ ν ∂ ν c i -ξ µ ∂ µ c i -χγ µ χA i µ -χγ µ χA i µ , (14)
E µ + E µ = -2χγ µ χ, G a + G a = 0. ( 15
)
Let us note that the introduction of an auxiliary real field Λ i besides the fields present in quantized N = 1 super Yang-Mills theory in four-dimensions, guarantees automatically the off-shell nilpotency of the {Q, Q} algebra and make easier then, as we will see in the next section, the gauge-fixing process.

III. QUANTUM ACTION

In the present section, we show how to construct in the context of our procedure a BRSTinvariant quantum action for N = 1 super Yang-Mills theory as the lowest component of a quantum superaction. To this purpose, we choose the following gauge-fixing superaction

S sgf = d 4 xL sgf , L sgf = (∂ 1 Φ 2 )(∂ µ Φ µ ) + (∂ µ Φ 2 )(∂ 1 Φ µ ) + (∂ 1 Φ 2 )(∂ 1 Φ 2 ). (16)
Let us recall that similar gauge-fixing superaction was used in Refs. [START_REF] Loumi-Fergane | [END_REF][20][21]. We note first that in the case of Yang-Mills theory the superaction involves a Lorentz gauge [22] given by

∂ µ Φ µ | = 0. (17)
In the case of super Yang-Mills theory we shall choose a supersymmetric gauge-fixing which is the extension of the Lorentz gauge. This gauge fixing can be obtained from ( 17) by using the following substitution

Φ µ → Φ µ = Φ µ + [∂ µ Φ a , Q a ]. (18)
Now, it is easy to see that the gauge-fixing superaction ( 16) can be put in the following form

S sgf = d 4 x[(∂ 1 Φ 2 )(∂ µ Φ µ ) + (∂ µ Φ 2 )(∂ 1 Φ µ )]. (19)
To determine the gauge-fixing action S gf as the lowest component of the gauge-fixing superaction S gf = S sgf | , we impose the following rules

T r(I m I n ) = δ m n , T r([Q a , Q b ]) = 2(γ µ ) ab ∂ µ , T r(P 2 ) = 0. ( 20
)
These rules permit us to compute the trace of each term in [START_REF] Loumi-Fergane | [END_REF] . Indeed, from (20) it is easy to put the gauge-fixing action Sgf in the form

Sgf = Ssgf | = d 4 x[B∂ µ A µ + 2b * j G(γ µ ∂ µ λ j ) +(∂ µ c){D µ c + ξ ν ∂ ν A µ + χγ µ λ} -2b * j (∂ µ χ)γ ν ∂ ν ∂ µ {f j ik λ i c k + ξ τ ∂ τ λ j - 1 2 χσ τ ν F j τ ν -χΛ j }]. (21)
On the other hand, the presence of the extrafield breaks the invariance of the classical action [START_REF] Henneaux | Quantization of Gauge Systems[END_REF]. In fact, the only terms which may contribute to the Q-variation of the classical action S 0 are those containing the extrafield Λ i . This follows from the fact that the BRST transformations up to terms Λ i represent the N = 1 super Yang-Mills transformations expressed à la BRST. A simple calculation with the help of the BRST transformations [START_REF] Bonora | [END_REF] leads to

QS 0 = χ a Λ i (γ µ ) ab (D µ λ b ) i . ( 22
)
Thus the classical action S 0 is not BRST-invariant, and in order to find the BRST-invariant extension S inv of the classical action, we shall add to S 0 a term S 0 so that

Q(S 0 + S 0 ) = 0. ( 23
)
We remark that S 0 is the part of the extended classical action related to the auxiliary field Λ i and is given by

S 0 = - 1 2 Λ i Λ i . (24)
Then, it is quite easy to show that Q(S 0 ) = -Q( S 0 ) by a direct calculation with the help of the transformations [START_REF] Bonora | [END_REF].

Having found the BRST-invariant extended action S inv we now write the full off-shell BRST-invariant quantum action S q by adding to the Q-invariant action, S inv = S 0 + S 0 , the Q-invariant gauge-fixing action S gf

S q = S 0 + S 0 + S gf . (25)
It is worth nothing that the quantum action (25) allows us to see that the auxiliary field Λ i does not propagate, as its equation of motion is a constraint

δS gf δΛ i = -Λ i + 2b * i ( χγ µ ∂ µ χ) = 0. ( 26
)
Thus the essential role of the nondynamical auxiliary field Λ i is to close the BRST and anti-BRST algebra off-shell.

The elimination of the auxiliary field Λ i by means of its equation of motion (26) leads to the same gauge-fixed theory with on-shell nilpotent BRST transformations obtained in the context of BV formalism [10] as well as in the framework of the superfibre bundle approach [START_REF] Loumi-Fergane | [END_REF].

Moreover, in our formalism we have also introduced an anti-BRST operator Q and it is important to realize that both the BRST symmetry and anti-BRST symmetry can be taken into account on an equal footing. To this end, we simply use the fact that there is a complete duality, with respect to the mirror symmetry of the ghost number, between the Q and Qtransformations. So, the Q-variation of the classical action S 0 is given by

QS 0 = χ a Λ i (γ µ ) ab (D µ λ b ) i . (27)
Using however the Q-transformations of the auxiliary field (see Eqs. [START_REF] Bonora | [END_REF] with the mirror symmetry), we obtain that the Q-invariant action S inv = S 0 + S 0 is also Q-invariant. Furthermore, the Q-gauge-fixing action can be also written as in Yang-Mills theories in Q-form.

Therefore the full off-shell BRST-invariant quantum action S q = S 0 + S 0 + S gf is also an off-shell anti-BRST-invariant quantum action.

IV. CONCLUSION

In the present paper we have developed a BRST superspace formalism in order to perform the quantization of the four dimensional N = 1 supersymmetric Yang-Mills theory as model where the classical gauge algebra is not closed. In this geometrical framework, our construction is entirely based on the possibility of introducing ab initio a set of fields through a super Yang-Mills connection. The latter represents the gauge fields and their associated ghost and antighost fields occurring in such theory, whereas the extrafield coming from the superconnection via the supersymmetric transformations is required to achieve the off-shell nilpotency of the BRST and anti-BRST operators. Let us note that for a local symmetry the minimal set of auxiliary fields is introduced by the supercurvature while for our case of global symmetry the auxiliary real field is introduced via the superYang-Mills connection.

Furthermore, we have performed a direct construction of the BRST invariant gauge fixed action for N = 1 , 4D supersymmetric Yang-Mills theory in analogy with what it is realized in BF theories [15] and simple supergravity [16]. The obtained quantum action allows us to see that the extrafield enjoy the auxiliary freedom. The elimination of this auxiliary field using its equation of motion permits us to recover the standard quantum action with onshell nilpotent BRST symmetry [10]. By using the mirror symmetry between the BRST and anti-BRST transformations, we can see that the BRST invariant extended classical action is also anti-BRST invariant. Therefore the full quantum action is BRST and anti-BRST invariant, since the gauge-fixing action can be written as in the Yang-Mills case in BRST as well as anti-BRST exact form, due to the off-shell nilpotency of the BRST-anti-BRST algebra.

Finally, we should mention that the BRST superspace formalism represents the natural arena where the fields and their off-shell nilpotent BRST and anti-BRST transformations for gauge theories can be found. This is not only the case of Yang-Mills type theories, arbitrary gauge theories may be also treated in this framework. Indeed, such formalism was applied to several interesting theories with local symmetry such non-Abelian BF theory [15] and simple supergravity [16]. In the present work this formalism has been applied successfully to the theory with a global supersymmetry. The off-shell nilpotency is naturally implemented through the introduction of auxiliary field required for the consistency of the BRST superspace geometry. Thus, it would be a very nice endeavor to use this basic idea to study the structure of auxiliary field in other gauge theories.

  ,4 are the Dirac matrices in the Weyl basis, b * i = b i for a = 1, 2 and b * i = -b i for a = 3, 4 giving the representation of the internal symmetry of Q a and [, ] the graded Lie bracket. Let us mention that the supersymmetric generators {Q a } are given in the Majorana representation

V. ACKNOWLEDGMENTS

A. Meziane would like to thank Prof. H. Nicolai, director at Max Planck Institute for Gravitational Physics, Albert Einstein Institute, Germany, for his warm hospitality and illuminating discussions. The DAAD is gratefully acknowledged too for his financial support.