Second-order phase-field formulations for anisotropic brittle fracture

We address brittle fracture in anisotropic materials featuring two-fold and four-fold symmetric fracture toughness. For these two classes, we develop two variational phase-field models based on the family of regularizations proposed by Focardi (Focardi, M. On the variational approximation of free-discontinuity problems in the vectorial case. Math. Models Methods App. Sci., 11:663{684, 2001), for which Gamma-convergence results hold. Since both models are of second order, as opposed to the previously available fourth-order models for four-fold symmetric fracture toughness, they do not require basis functions of C1-continuity nor mixed variational principles for finite element discretization. For the four-fold symmetric formulation we show that the standard quadratic degradation function is unsuitable and devise a procedure to derive a suitable one. The performance of the new models is assessed via several numerical examples that simulate anisotropic fracture under anti-plane shear loading. For both formulations at fixed displacements (i.e. within an alternate minimization procedure), we also provide some existence and uniqueness results for the phase-field solution.


Introduction
Predicting the fracture behavior in brittle materials with anisotropic fracture toughness is still a challenge. The extension to the anisotropic case of selection criteria for the crack direction which have proved successful in fracture mechanics of isotropic materials, such as the principle of local symmetry [22,12] and the maximum energy release rate criterion [28], is by no means trivial [34]. The generalized energy release rate criterion [37,24,34], where the crack propagation direction is chosen by maximizing the ratio of the energy release rate to the fracture toughness with respect to the propagation angle, was shown to correlate successfully with experiments on materials exhibiting both weak [29] and strong [39] anisotropy in the fracture toughness, see also the detailed review and discussion in [33].
The variational phase-field approach to fracture, pioneered by Bourdin et al. [5] and first proposed as the regularization of Francfort and Marigo's variational fracture formulation [17], was recently extended to materials with anisotropic fracture toughness. Hakim and Karma [26,27] proposed a phase-field model for materials featuring two-fold weak anisotropy and elucidated several features of crack propagation in these materials. This model was extended to the three-dimensional setting including geometric nonlinearity and nonlinear elasticity [11] and to the consideration of multiple cleavage planes using multiple phase-field variables [36]. With the purpose of tackling phenomena such as sawtooth crack patterns and forbidden crack directions, which are directly related to four-fold strongly anisotropic surface energies, Li et al. [32] formulated a fourth-order anisotropic variational fracture model inspired by phase-field models of crystal growth. Similarly, Teichtmeister et al. [40] formulated second-and fourth-order phase-field models based on structural tensors for twofold and four-fold symmetry, respectively, and extended the framework to large deformations. Li and Maurini [33] improved and simplified the phase-field model in [32] in order to enable an analytical solution for the optimal crack profile, while keeping the fundamental properties of the phase-field approximation. Once again, four-fold symmetry was addressed with a fourth-order model.
In this paper, we address brittle fracture in anisotropic materials featuring two-fold and four-fold symmetric fracture toughness and develop two variational phase-field models based on the family of regularizations proposed by Focardi [16], who also proved their Γ-convergence. The main novelty lies in the model for four-fold symmetric anisotropy: as opposed to all previously available models for the same type of anisotropy, our proposed model is of second order, hence it does not require higher-continuity basis functions nor mixed variational principles for finite element discretization. For the four-fold symmetric formulation we show that the standard quadratic degradation function is unsuitable and propose a procedure to derive a suitable one. The performance of the new models is assessed via several numerical examples that simulate anisotropic fracture under anti-plane shear loading. For both formulations at fixed displacements, we also provide some existence and uniqueness results for the phase-field solution.
The remainder of this manuscript is organized as follows. Section 2 introduces basic concepts on anisotropy of fracture properties that serve as basis for the following investigation. In Section 3, we recall the phase-field formulation of (isotropic) brittle fracture based on the Ambrosio-Tortorelli regularization and introduce the main results in [16], which are our point of departure for the new developments in the subsequent sections. We show that the models based on the Focardi regularization combined with the Euclidean norm, which we term isotropic Foc-p models, can be viewed as a family of models which includes the wellknown AT-2 model as a special case (corresponding to the isotropic Foc-2 model). In Section 4, we focus on the isotropic Foc-4 model. We show that the standard quadratic degradation function is unsuitable for this model, we derive a new suitable one, and we test the ensuing formulation on an anti-plane shear numerical setup. In Section 5, we develop the anisotropic counterparts of the isotropic Foc-2 and Foc-4 formulations, which model an anisotropic fracture behavior with two-and four-fold symmetric fracture toughness, respectively. Section 6 illustrates the numerical tests on the two anisotropic formulations for different anti-plane shear setups. We draw the main conclusions in Section 7. Appendix A reports some existence and uniqueness results on both anisotropic formulations.
2 Basic concepts on anisotropic fracture toughness :probsetting Materials with anisotropic fracture behavior are characterized by an orientation-dependent fracture toughness or critical energy release rate G c . One of the examples is sketched in Figure 1a, where it is assumed that there exists a plane of anisotropy A ⊂ R 2 within the body Ω ⊂ R 3 such that for any point x ∈ A and with the polar coordinate system within A associated to x, G c is viewed as a function of the polar angle θ ∈ [0, 2π), that is, G c = G c (θ). This is in contrast to the simple isotropic case when G c ≡ const. In this paper, we consider the two-dimensional setting by taking Ω ⊂ R 2 , which implies that the anisotropy plane A and Ω coincide. For modeling purposes, it is furthermore assumed that where G 0 is a positive constant and γ(θ) is a periodic trigonometric function. A suitable form of this function reads The parameter ω defines the direction corresponding to the largest fracture toughness and hence it is referred to as the strongest or principal material direction. The parameter τ defines the magnitude of the largest fracture toughness and is called the anisotropy strength. A useful induced quantity for characterization of the anisotropy strength is the ratio 1+τ 1−τ ∈ [1, ∞). Finally, k is the periodicity parameter such that the function γ k is 2π k -periodic. The above gives rise to the notion of a k-fold symmetric fracture toughness G c .
In the present work, we are interested in the cases k = 2 and 4 in (2), i.e. in the so-called two-and four-fold anisotropy cases. We give the corresponding explicit expressions of γ k for future reference: with τ ∈ [0, 1), ω ∈ [0, π), and γ 4 (θ) = 1+τ cos(4(θ − ω)) = 1+τ cos(4ω)(cos 4 (θ) − 6 sin 2 (θ) cos 2 (θ) + sin 4 (θ)) . The polar plots of both functions are continuous closed curves surrounding the origin of the polar coordinate system. According to standard terminology, the anisotropy is said to be weak if γ (θ) is convex and strong if γ (θ) is non-convex. The convexity of the function γ (θ) is known to be equivalent to the convexity of the region enclosed by the polar plot of γ −1 (θ), which is guaranteed if With the expressions of γ 2 (θ) and γ 4 (θ) in (3) and (4), it is straightforward to verify that regardless of the value of ω. The polar plots of γ 2 (θ) and γ −1 2 (θ), as well as γ 4 (θ) and γ −1 4 (θ) for a fixed ω and different magnitudes of τ ∈ [0, 1) are depicted in Figure 1b. With the increase of τ , the elliptic shape of γ 2 and γ −1 2 turns into a peanut-like shape, whereas for γ 4 and γ −1 4 one has a transformation from super-ellipse to 'butterfly'.  In the following, when considering the cases of isotropic and anisotropic fracture toughness we refer to them as 'isotropic fracture' and 'anisotropic fracture', respectively. In the former case, we assume G c (θ) ≡ const such that, owing to (1) where γ(θ) ≡ 1, it is simply G c (θ) = G 0 . In the latter case, the term 'anisotropic fracture' is accompanied by the specification 'with k-fold symmetric G c '.
3 Variational formulation of brittle fracture and phase-field regularization sec:varform In this section, we briefly recall the phase-field formulation of (isotropic) brittle fracture from [5,6,7,8], based on the so-called Ambrosio-Tortorelli regularization [2,1]. This is followed by the presentation of the main results by Focardi [16], which are our point of departure for the new developments in the subsequent sections. Finally, we briefly outline the incremental minimization problem which is solved for computing a quasi-static phase-field evolution.
Let Ω ⊂ R d , d = 2 or 3 be an open and bounded domain representing the configuration of a d-dimensional body, and let Γ D,0 and Γ D,1 be the (non-overlapping) portions of the boundary ∂Ω of Ω which are fixed and subjected to the imposed displacementū, respectively ( Figure  2a). Let also Γ c ⊂ Ω represent the crack set, that is, the set of points where the displacement field u is discontinuous.

AT_reg
The variational approach to (isotropic) brittle fracture by Francfort and Marigo [17] relies on the energy functional VAF and the related incremental minimization problem. In (8), Ψ is the elastic energy density, a function of the infinitesimal strain ε(u) = 1 2 (∇u + (∇u) T ), where u : Ω → R d is the displacement field and ∇ is the gradient operator, and H d−1 is the Hausdorff measure. The terms E el. and E S respectively represent the elastic energy stored in the cracked body and the surface energy associated to the crack. Upon integration, the latter reads E S (Γ c ) = G 0 H d−1 (Γ c ) where, in simple terms, H 1 (Γ c ) and H 2 (Γ c ) are the length and the surface area of Γ c when d = 2 and 3, respectively. Thus the expression of the surface energy associated to the crack is consistent with Griffith's theory [23] with a fracture toughness G 0 .
The regularization of (8)à la Ambrosio-Tortorelli [2] developed by Bourdin-Francfort-Marigo [5,6,7,8], which forms the basis for a variety of (isotropic) fracture phase-field formulations, reads as follows: with α : Ω → [0, 1] as the phase-field variable, which takes the value 1 on Γ c , decays smoothly to 0 in a subset of Ω\Γ c and vanishes in the rest of the domain, as sketched in Figure 2b. With this definition, the limits α = 1 and α = 0 represent the fully broken and the intact (undamaged) material phases, respectively, whereas the intermediate range α ∈ (0, 1) mimics the transition zone between them. The function g, often denoted as the degradation or the coupling function, is responsible for the material stiffness degradation, the function w, often termed the local dissipation function [38], defines the (decaying) profile of α, whereas the parameter 0 < diam(Ω) controls the thickness of the localization zone of α, i.e. of the transition zone between the two material states.
The specific choice of the functions g and w in (9) establishes the rigorous link between (8) and (9) when → 0 via the notion of Γ-convergence, see e.g. [9,10,1], also giving a meaning to the induced constant c w . Thus, g is a continuous monotonic function that fulfils the properties: g(0) = 1, g(1) = 0, g (1) = 0 and g (α) < 0 for α ∈ [0, 1), see e.g. [38]. The quadratic polynomial g(α) := (1 − α) 2 is the simplest choice. The function w is continuous and monotonic such that w(0) = 0, w(1) = 1 and w (α) ≥ 0 for α ∈ [0, 1]. The constant c w := 4 1 0 w(t) dt is a normalization constant in the sense of Γ-convergence. The two suitable candidates for w which are widely adopted read w(α) ∈ {α, α 2 }, such that c w ∈ { 8 3 , 2}, respectively. The two instances of formulation (9) containing the aforementioned choices for g and w are typically termed AT-1 and AT-2 models, see Table 1. AT stands for Ambrosio-Tortorelli type of regularization, see [2,1]. Table 1: Ingredients of formulation (9). g w c w Name for (9) The main difference between the two models is that AT-1 leads to the existence of an elastic stage before the onset of fracture, whereas using AT-2 the phase-field variable starts to evolve as soon as the material is loaded, see e.g. [38,3,35] for a more detailed explanation, as well as the discussion in Section 4.1.

Foc_reg
In Focardi [16], the Ambrosio-Tortorelli regularization of the free-discontinuity problem is generalized so as to include an orientation-dependent surface energy, and the related Γconvergence result is obtained. Using the notation in (8) and (9), the unregularized energy Table 2: Ingredients of formulation (11) assuming ϕ(a, b) = (a 2 + b 2 ) To be specified (Section 4) Isotropic Foc-4 model Foc functional considered by Focardi and its regularized (phase-field) counterpart read RegVAF_f respectively. In the above, n is the normal unit vector to Γ c , ϕ : The coupling function g is supposed to fulfill the standard properties listed in Section 3.1. Theorem 3.1 in [16] establishes the Γ-convergence result between E(u, α) and E(u, Γ c ) under no specific assumptions for the norm ϕ.
For any choice of ϕ, one may say that (11) provides a p-parametric family of regularizations to the formulation in (10). We denote this family as the Foc-p models.
Suppose now that ϕ is a Euclidean norm. In this case, since n is a unit vector, we simply arrive at ϕ(n) ≡ 1 in (10) -this then becomes identical to (8) -and we also obtain ϕ p (∇α) = |∇α| p in (11). In the following, we term the p-parametric family of formulations in (11) with ϕ given by the Euclidean norm the isotropic Foc-p models (when the reference to isotropy is clear from the context, the term isotropic is omitted).
The isotropic Foc-p model with p > 1 can be viewed as an alternative to the Ambrosio-Tortorelli approximations AT-1 and AT-2, and, in particular, a generalization of the AT-2 case. In Table 2 the isotropic Foc-2 and Foc-4 models are reported. One immediately notices that the Foc-2 formulation is identical to the AT-2 model, provided that the degradation function in (11) is taken as g(α) = (1 − α) 2 . The Foc-4 formulation turns out to be a more interesting case and is thoroughly investigated in Section 4.
However, the generality of (10) and (11) is significantly broader, since ϕ does not have to be a Euclidean norm. In Section 5, we explain how the norm ϕ in (10) -hence in (11) -as well as the parameter p in (11) can be chosen based on the anisotropy functions γ k , k ∈ {2, 4} from (2), and thus used to model anisotropic fracture (we will find out that the most natural choice is p = k). This results in what we term the anisotropic Foc-p models. Table 3 summarizes the terminology introduced so far.

Incremental variational problem
We consider a quasi-static loading process in which the monotonically increasing imposed displacementū n , with n = 1, 2, ... denoting the discrete pseudo-time step parameter, is  (10) and (11) Name for (11) Aims at phase-field modeling of Not specified Foc-p models -

Euclidean
Isotropic Foc-p models Isotropic fracture γ p -induced Anisotropic Foc-p models Anisotropic fracture with p-fold symmetric G c Mnem prescribed on Γ D,1 . It is assumed that, during this process, the crack surface Γ c evolves incrementally and in an irreversible manner, that is, Γ c,n ⊇ Γ c,n−1 .
With the phase-field energy functional E defined by (9) or (11), the state of the system at a given loading step n ≥ 1 is given by Here is the space of kinematically admissible displacements at step n, with W 1,2 (Ω; R d ) as the usual Sobolev space of functions with values in R d , and is the admissible space for α at step n, with W 1,p (Ω; R) as the usual Sobolev space of functions with values in R and α n−1 known from the previous step. The condition α ≥ α n−1 in Ω is used to enforce the irreversibility of the crack phase-field evolution. It is the backward difference quotient form ofα ≥ 0 in Ω. Also, in D α n−1 , we take p = 2 for the AT-1 and AT-2 models in (9), and p > 1 for the Foc-p model in (11). The necessary optimality conditions for (u, α) ∈ V un × D α n−1 at every loading step n ≥ 1, written in partitioned form, read: where E u and E α denote the directional derivatives of E with respect to u and α, respectively. The displacement test space in (13) is Conditions (13) characterize, in general, a local minimum of E.

Rem_Irr
In what follows, we incorporate the irreversibility constraint α ≥ α n−1 via simple penalization 1 , thus arriving at RegVAF_pen where E is given by (9) or (11), a − := min(0, a) andλ ∈ R + is the penalty parameter. In [20], an analytical procedure for the 'reasonable' choice of a lower bound forλ that guarantees a sufficiently accurate enforcement of the constraint is devised. The optimalλ is given bŷ lambda_opt where 0 < TOL ir 1 is a user-prescribed irreversibility tolerance threshold (in [20], a practical value for TOL ir is suggested as 0.01).
With F defined by (14), the sought solution at a given loading step n ≥ 1 is given by The optimality condition in this case simplifies to where F u = E u , and Notice that in (16) and (17) we again take p = 2 for the AT-1 and AT-2 models, and p > 1 for the Foc-p model.

Res1
In Section 3.2 we have seen that the isotropic Foc-p models with p > 1 can be viewed as a family of models which includes AT-2 as a special case. In view of the future extension to anisotropic fracture (to be addressed in Section 5), and in particular to the two-and four-fold anisotropy cases, we are especially interested in the isotropic Foc-2 and Foc-4 models. Since the former (with the standard quadratic degradation function) is identical to the AT-2 model, whose properties are already well known, in this section we concentrate on the isotropic Foc-4 model (which we simply denote as Foc-4 model throughout this section). The corresponding energy functional is obtained from (11) by setting ϕ(a, b) := (a 2 + b 2 ) 1 2 and p = 4, which yields where b w = 2 − 4 3 . As follows, we first carry out the analysis of a one-dimensional bar under homogeneous conditions. Through this analysis we show that the standard quadratic degradation function is not suitable for the Foc-4 formulation, and we then derive analytically a new suitable expression for g. Subsequently, we study the one-dimensional localization behavior and construct the so-called optimal profile of α induced by the fracture energy term E S (α) in (18). We compare it with the optimal profiles of the AT-1 and AT-2 models. The knowledge of this profile also allows us to estimate the magnitude of the penalty parameterλ in (14) to be used for the enforcement of irreversibility in computations. Finally, the new formulation obtained combining the Foc-4 model and the new degradation function is tested on a simple anti-plane shear numerical setup.

Res1-g
The main tool for assessing the fundamental properties of a fracture phase-field formulation -in particular, the impact of ingredients such as g -is the analysis of a one-dimensional elastic bar under tension [38,35,4,31]. In the following, we recall the corresponding results for the AT-1 and AT-2 models, and apply the same analysis to understand the role of the degradation function and to devise a suitable one for the Foc-4 model.
Following [38], we consider a one-dimensional elastic bar of length L, clamped at one end (u(0) = 0) and subjected to an imposed displacementū at the other end (u(L) =ū). The elastic energy density function reads Ψ( ) = 1 2 E 2 with E as the Young's modulus. In this pure tensile loading case, the phase-field irreversibility constraint is not necessary, hence we can setλ = 0. The strong form of the equilibrium equations in (0, L) reads (regardless of the model), whereas the damage criterion reads for the AT-1 and AT-2 models, and for the Foc-4 model. Here, we use the re-scaled quantities¯ = / G 0 E andσ = σ/ G 0 E . Let us focus on the study of the homogeneous solution, where the damage and the stress and strain fields are uniform along the bar [38,35]. In this case, due to the disappearance of the spatial derivatives of α, equations (19)-(21) provide the explicit¯ -α andσ-α relations: • for the AT-1 model,¯ .
In the first two cases, we already substituted for g the standard quadratic expression, whereas we are considering a generic g for the third case. Figure 3 depicts the¯ -α andσ-α plots for the AT-1 and AT-2 models, as well as for the Foc-4 model with g(α) = (1 − α) t , t ∈ {2, 3, 4}.  The curves obtained with the AT-1 and AT-2 models are well known. From them it can be inferred that the AT-1 model features an initial linearly-elastic phase of behavior, after which damage starts to occur. The peak value of the stress is reached for a value of α given by α cr = 0. Conversely, the AT-2 model has no purely elastic phase, as the evolution of α starts from the onset of loading and the peak value of the stress is reached for α cr = 0.25.
The curves corresponding to the Foc-4 model are qualitatively similar to those of the AT-2 model. However, the non-linearity prior to the peak stress is even more pronounced and the value of α cr is even larger. Indeed, for g(α) = (1 − α) 2 it is α cr = 0.5. If the polynomial degree of the degradation function is increased, α cr is reduced, however α reaches 1 at increasingly large values of¯ .
We now aim at devising a degradation function for the Foc-4 model that, in addition to fulfilling the basic properties of degradation functions (see Section 3.1), leads to a¯ -α behavior similar to the one of the AT-1 model. We thus enforce the following property 2 where s ∈ R and m ∈ N, m ≥ 1. Solution of the differential equation in (25), with the boundary conditions g(0) = 1 and g(1) = 0, reads Secondly, forσ given by (24) with g(α) as in (26), we seek m ≥ 1 to yield dσ dα ≤ 0 in α ∈ [0, 1]. Figure 4 illustrates dσ dα as a function of (α, m). It can be clearly seen that the desired property for this quantity is fulfilled for m = 4, yielding (27) g_8 Figure 5 shows the¯ -α,σ-α andσ-¯ relations for the Foc-4 model with the degradation function (27), also in comparison with the AT-1 and AT-2 cases. It is clear that the constructed eighth-order polynomial function in (27) leads to the desired behavior.

Res1-profile
If Γ c ⊂ Ω is a fully developed crack, the phase-field profile α : Ω → [0, 1] that represents the regularized counterpart of Γ c can be obtained by minimizing the fracture surface energy functional E S on a suitable admissible set, namely Optim see Figure 6a for an illustrative sketch in two dimensions. Using a slicing argument, one may resort to a one-dimensional setting in (28): we assume Ω := (−L, L) with L such that Γ c is represented by the point t = 0, and we also denote α := α(t). In this setting, the minimizer α can be constructed explicitly.
The strong formulation of the boundary value problem to compute α for the AT-1 and AT-2 models is given by Note that in (29) the inequalities hold for AT-1 and the corresponding equalities for AT-2. The solutions representing the optimal phase-field profiles read for AT-1, and for AT-2, see e.g. [20] for details.
For the Foc-4 model, the optimal profile must be computed from and is straightforwardly obtained as 3 The plots of α in (30), (31) and (33) are shown in Figure 6b. It is evident that the AT-2 and the Foc-4 models have very similar optimal profiles. A practical outcome of this observation is that, when using the penalty method to enforce irreversibility, we can choose the penalty parameterλ for the Foc-4 formulation using the same estimate valid for the AT-2 formulation, see eq. (15).
The applied incremental displacement is given byū n = n∆ū, n = 1, ..., 3 2∆ū , with ∆ū = 0.1 as the loading increment. We employ the numerical package FreeFem++ [25]. Both the displacement field u and the crack phase-field α are approximated using P 1 -triangles. The finite element mesh is a pre-adapted one, which is refined in the region of Ω where crack propagation is expected. The minimum and maximum characteristic mesh sizes are given by (h min , h max ) := ( 1 5 , ), where h min and h max stand for the mesh size inside and outside of the refined region, respectively. Figure 8a-b visualizes the α fields computed with the AT-2 and Foc-4 models with quadratic degradation function at the same loading step after the peak load. It is evident that, in both cases, α does not vanish outside of the localization zone, which is a direct consequence of the accumulation of damage since the onset of loading discussed in Section 4.1. For the given setup at the chosen loading step, it is min(α) ≈ 0.05 for AT-2 and min(α) ≈ 0.25 for Foc-4. Interestingly, both values are significantly lower than the value of α cr obtained for the one-dimensional homogenous responses of the respective models, see Section 4.1. Finally, Figure 8c shows results for the Foc-4 model with the degradation function in (27), where no spurious damage accumulation outside of the localization region is obtained.  In this section, we develop the anisotropic counterparts of the isotropic Foc-2 and Foc-4 formulations. To this end, we suitably choose the norm ϕ in (10) -and, as a result, in (11) -based on the knowledge of the anisotropy functions γ k (θ), k ∈ {2, 4} in (3) and (4). In particular, γ 2 is used to induce the norm in the Foc-2 formulation, whereas γ 4 is employed for defining ϕ in the Foc-4 formulation. We term the resulting formulations anisotropic Foc-2 and Foc-4 models, see Table 3.

Derivation of the anisotropic models
Since ϕ in (10) is assumed to be a function of the unit normal n = (n x , n y ), whereas the anisotropy functions γ k take as argument the polar angle θ, the first straightforward step is to express θ via the components of n. Let Γ represent a pre-existing crack, point x ∈ Γ be the coordinate of its tip, and Γ ext be a possible extension of Γ starting from x. Also, let us define the Cartesian and polar coordinate systems associated to x, as well as the unit normal n to Γ ext at x. The angle θ is defined as the angle between the positive direction of the x-axis and the crack extension Γ ext , see Figure 9. With this definition, θ is related to the components of n by n x = cos(θ + π 2 ) = − sin(θ) and n y = sin(θ + π 2 ) = cos(θ). Substituting in (3) and (4), we obtain γ 2 (n) = n 2 x + n 2 y + τ cos(2ω)(n 2 y − n 2 x ) − 2 sin(2ω)n x n y , (34) gamma2-n and γ 4 (n) = (n 2 x + n 2 y ) 2 + τ cos(4ω)(n 4 y − 6n 2 x n 2 y + n 4 x ) − 4 sin(4ω)n x n y (n 2 y − n 2 x ) .
where we also exploited the following representations of the unity

unity
We can now introduce the mappings ϕ : R 2 → [0, +∞) for two-and four-fold anisotropy, which respectively read: We call ϕ in (37) and (38) γ 2 -and γ 4 -induced norms, respectively, to be used in the Focardi discrete variational formulation in (10). Noting now that the above norms ϕ to be used in (10) have the form ϕ = f 1 k , k ∈ {2, 4}, and that the consequent terms ϕ p , p > 1 to appear in (11) read ϕ p = f p k , the most natural choice is to set p ! = k ∈ {2, 4}. This yields as the ingredient of (11) when p = 2 and ϕ is the γ 2 -induced norm, and as the ingredient of (11) when p = 4 and ϕ is the γ 4 -induced norm. We denoted (α ,x , α ,y ) as the components of ∇α. The norm in (39) can be also written as with i, j = 1, 2. Here is a symmetric positive definite second-order tensor, and I is the two-dimensional second-order unit tensor. The norm in (40) can be written as ϕ 4 (∇α) = B∇α ⊗ ∇α · ∇α ⊗ ∇α = B ijkl α ,i α ,j α ,k α ,l (43) phi4_tensor with i, j, k, l = 1, 2. Here is a positive definite fourth-order tensor exhibiting minor and major symmetries, with components   (37) and (38), respectively -is the only appropriate one in this case, as it preserves an identical dimension for all terms entering the corresponding ϕ.

Remark 2
The choice of the exponents 1 2 and 1 4 in (37) and (38) is a natural one, since the limiting case of τ = 0, with (a, b) = 0, turns the corresponding ϕ into the Euclidean norm. (37) and (38) indeed define two norms, since owing to

Remark 3 The mappings in ϕ in
the following norm-equivalence result holds: for all τ ∈ [0, 1) and ω belonging to the appropriate range.
The questions of existence and uniqueness of the solution to the minimization problem for E in (51) and (52) is a delicate task. Indeed, it is very well known that for the isotropic case, min (u,α) E may admit infinitely many solutions, or none. In the light of the so-called alternate minimization scheme, see e.g. [5,6,7,8], applied to resolve min (u,α) E numerically, one typically restricts the discussion to studying the well-posedness of the problems min u E when α is fixed, and min α E when u is fixed. The former case is simple: the problem is wellposed due to the convexity of E with respect to u. Well-posedness of the latter minimization problem with E given by (51) and (52) is discussed in Appendices A.1 and A.2, respectively. In the appendices, this is formulated as an τ -parametric result.

Relation to existing phase-field models for anisotropic fracture
It is straightforward to recognize that the surface energy integral in (51) is a specific instance of the class of integrals proposed in the literature for two-fold anisotropy [26,27,33] where w(α) and c w are specialized to the AT-2 model. It is more interesting to compare the surface energy integral in (52) with the integral proposed in [33] for four-fold anisotropy (a simplified version of the model in [32]): Here the fourth-order tensor B is a positive definite fourth-order tensor exhibiting minor and major symmetries, just like the B tensor in (43). However, our new formulation has the advantage of involving only first-order derivatives of α, so that a standard C 0 finite element approximation can be adopted, whereas the second-order derivatives of α contained in (54) call for C 1 -continuous finite element approximations or mixed approaches.

Numerical solution
Let the energy functional F be as in (14) with E represented by (18), (51) and (52). The staggered iterative process for solving F(u, α) → min at any fixed loading step n ≥ 0 consists in the following three steps: with the known initial guess (u 0 , α 0 ), for any staggered iteration k ≥ 1, 3. for the computed pair (u k , α k ) check As shown in Appendix A.2, both the isotropic and anisotropic Foc-4 formulations in (18) and (52), respectively, are non-convex with respect to the crack phase-field variable α. The implication is that solving F α (u k−1 , α; β) = 0 with standard Newton-Raphson iterative procedures may lead to convergence difficulties. In the context of phase-field modeling, a sophisticated line search algorithm was developed in [19]. In this paper, instead of applying the same algorithm for the case at hand, we investigate the use of the trust-region (TR) method [13], first adopted in [30] in combination with phase-field simulations.
In the version we use, the TR method relies on the second-order Taylor expansion of F around a given fixed α: with F αα as the second-order directional derivative of F with respect to α. The iterative procedure of obtaining α k is as follows: given α (i) , i ≥ 0 (as the initial guess, one may take α (0) := α k−1 ), the minimization problem to be solved at every i reads Here, R i > 0 is the TR radius, considered as an iteration-dependent quantity. Notice also that recent m must not be confused with m from section 4.1. For the computed z, we obtain the trial step solution α (i+1) = α (i) + z and the corresponding TR ratio If 0 < η 1 < ρ i < η 2 < 1, the solution α (i+1) is considered acceptable and the TR radius R i does not need to be modified. If ρ i ≤ η 1 (it can happen to be even negative), the solution is considered unacceptable, R i is reduced by factor γ 1 < 1 and z has to be recomputed. Finally, when ρ i ≥ η 2 , α (i+1) is acceptable and R i can be enlarged by factor γ 2 . Typical values for the parameters are η 1 = 1 4 , η 1 = 3 4 and γ 1 = 1 2 , γ 2 ∈ {2, 3}. The iterative solution process in i is terminated once a certain criterion is achieved, and we set the last α (i) to be α k .
Further relevant details are highlighted as follows: • The minimization in (55) is carried out under the point-wise constraint −R i ≤ z ≤ R i in Ω. This is the strongest constraint of the type ||z|| L m (Ω) ≤ R i , m ≥ 1, but is also the simplest one in terms of implementation. We incorporate it via penalization by minimizing the modified functional: MiPen with λ 1. In our numerical experiments, we typically take R 0 = 0.01 (this proves reasonable since the crack phase-field variable is bounded α ∈ [0, 1], whereas z can be viewed as a 'correction' and hence cannot be 'large'). Also, in all following examples we set the penalty parameter λ in (57) as 10 4 .
• To compute z = arg min{ M i (z) : z ∈ H 1 (Ω)} at every fixed i ≥ 0, the optimality Due to the Macaulay brackets in the penalty term, the problem is non-linear in z and can be solved e.g. with the Newton-Raphson method.
One important remark here is that the initial guess z 0 for the Newton-Raphson update z 0 + ∆z should not be 0 in order not to solve the original badly-behaved problem F α (u k−1 , ∆z; Y ) = 0.
• In (56), M i is substituted by M i .
• Our termination criterion for the acceptable trial step solution α (i+1) reads as follows. First, let us notice that checking whether the residual F α (u k−1 , α (i+1) ; β) is sufficiently small is an inappropriate criterion due to non-convexity of F in α: even if F(u k−1 , α (i+1) ) < F(u k−1 , α (i) ) holds true, as required, it can very well happen -and we indeed observed this in the computations -that F α (u k−1 , α (i+1) ; β) F α (u k−1 , α (i) ; β). As a result, we opt for monitoring and checking the relative error in F and stop iterating in i when the condition ≤ TolPF is fulfilled. In our experiments, we use TolPF = 10 −4 .

Numerical experiments
NumTest As follows, we illustrate some numerical experiments on the presented anisotropic Foc-2 and Foc-4 formulations. Given that the former belongs to an already known family of models for two-fold anisotropy, we mainly concentrate on the latter, which for the first time describes four-fold anisotropy with a second-order model. We restrict ourselves to an anti-plane shear setup. In addition to being computationally less expensive due to a scalar unknown displacement field, this setup offers the advantage that the anisotropy of the elastic stiffness properties does not play a role, at least as long as we consider anisotropy classes where normal and shear stresses remain uncoupled (e.g. orthorombic or cubic).
The applied incremental displacement is given byū n = n∆ū with ∆ū = 0.1 and n = 1, ..., 3 2∆ū . Both u and α are approximated using P 1 -triangles. The finite element mesh is identical for all computations, which only differ by the value of τ . The mesh is pre-refined in the region where crack propagation is expected and the characteristic mesh sizes in the refined and in the coarser regions are (h min , h max ) := ( 1 5 , ). The results are presented in Figure 10, where for comparison we also show results obtained with the isotropic Foc-2 (≡AT-2) formulation. It can be observed that, for a non-zero anisotropy strength τ , the crack path deviates from the vertical direction obtained in the isotropic case. As expected, for all considered values of τ , the crack tends to propagate along the weakest material direction (in this case, at the angle − π 4 ), and the increase of τ makes this effect more pronounced: already for τ = 0.5 (such that 1+τ 1−τ = 3) the angle of propagation is very close to − π 4 at the initial loading stage and then changes due to boundary effects, whereas for the largest τ = 0.8 (yielding 1+τ 1−τ = 9) the propagation angle aligned with the weak direction − π 4 persists until complete failure. Figure 10: Phase-field contour plots at the final loading step for the setup in Figure 7. The anisotropy function is γ 2 (θ) with the specified values of (ω, τ ), as depicted in the lower inserts. Fig:s0_Foc2 6.2 Numerical experiments for the anisotropic Foc-4 model

Res2-anis4
In this section we perform numerical tests for the anisotropic Foc-4 formulation in (52), using the anti-plane shear loading setups in Figures 11 (example 1       In example 1, the domain consists of a central strip of material featuring four-fold anisotropic fracture toughness, and two side strips of isotropic material with G c (θ) ≡ 100G 0 . For simplicity, the shear modulus µ is identical for both materials. For this example, the strongest material direction ω and the anisotropy strength parameter τ are fixed to 0 and 0.5, respectively. In examples 2 and 3, the domain consists entirely of anisotropic material, and for the fixed parameter τ = 0.5 we compare two values of ω ∈ {0, π 4 }.

Example 1
In this example, the idea of the bi-material arrangement with strong direction ω = 0 is to induce a crack evolution along the weakest material directions within the middle strip, accompanied by crack reflections at the boundaries with the two side strips (where the crack cannot propagate due to the significantly larger fracture toughness).
The results are presented in Figure 14, where the crack phase-field is computed on meshes with different element sizes. The first observation is that the actual crack propagation pattern only partially corresponds to the expectations. Let us start by looking at results for (h min , h max ) := ( 1 4 , ). The crack does start propagating along the weak material direction − π 4 . However, upon hitting the interface with the side strip of very tough material, it is not immediately deviated to the symmetric weak material direction π 4 , but rather proceeds for a certain length along the bimaterial interface and then deviates to π 4 . When the crack hits the other bimaterial interface, once again no immediate reflection takes place, and the crack proceeds along the bimaterial interface until reaching the end of the specimen.  A second observation stems from the comparison between the results obtained with the three meshes. It appears that the crack patterns are qualitatively similar. Also, it is evident that the results obtained for the last two meshes are roughly symmetric with respect to those of the first mesh. Due to symmetry, for each crack evolution two symmetric patterns are fully equivalent, and which one is obtained numerically is determined by the mesh, as well as by parameters related to the solution algorithm (including e.g. loading increments, thresholds and tolerances), and by round-off errors. Note that multiple solutions are possible already in isotropic fracture, see [21], but the range of possibilities appears to be even wider in the anisotropic case.
A third observation is the following. As mentioned above, we expected the crack to be reflected immediately upon hitting the bimaterial interface, and not to follow it as seen in . We illustrate how the crack does not reflect from the interface along one of the weakest material directions, as expected, but propagates along the strongest material direction in a micro zig-zagging mode.  Figure 14. Indeed, the bimaterial interface is aligned with the strongest material direction. Interestingly, even though the crack macroscopically propagates in this 'prohibited' direction, at a closer look it is seen to consist of many small segments following the two weak material directions. We refer to this phenomenon as 'microscopic (or micro) zig-zagging' and evidence it in Figure 15. Here the term 'microscopic' is intended to refer to a size scale of the order of the size of a finite element, as opposed to the macroscopic scale which relates to the characteristic size of the domain under investigation.
In order to observe whether the number of reflections at the bimaterial interface is influenced by the length-to-width ratio of the middle strip, we numerically test two additional specimens where the middle strip is either narrower or longer than in the specimen of Figure  14, see Figure 16. Indeed, when the length-to-width ratio of the strip is increased by either decreasing its width or increasing its length, the number of reflections increases by one. Also in these patterns, the crack partially follows the bimaterial interface between subsequent reflections and along the final part of the propagation, which locally induces once again a micro zig-zagging phenomenon.

Example 2
In this example, due to the presence of two initial cracks, we expect both cracks to propagate and to join in the middle of the specimen, following a pattern depending on the choice of the strongest material directions. Results are shown in Figure 17 for ω = 0 and ω = π 4 . In the first case, the simplest expected crack pattern would consist in the propagation of both initial cracks along + π 4 (for the left crack) and − π 4 (for the right crack), until both cracks join on the middle axis of the specimen. The obtained results are similar to this expectation, with the exception of a final segment followed by both cracks in the horizontal (i.e. the strongest) material direction, and again featuring a micro zig-zag pattern. In the second case, the two cracks propagate following first the vertical and then approximately the horizontal direction, i.e. both weakest directions. However, also in this case some local deviations are observed, which break the symmetry of the results and introduce again micro zig-zagging patterns. These local deviations and the loss of symmetry intuitively suggest that many alternative crack patterns could be obtained through small perturbations of the system, see [21].    Like example 2, this example features two initial cracks, which however are situated on opposite sides of the specimen. Similarly as in the previous two examples, we expect the two cracks to propagate along the weakest material directions and to merge at some point in the middle of the specimen. Once again, we examine both ω = 0 and ω = π 4 . Results are illustrated in Figure 18.
For ω = 0, the two initial cracks propagate approximately along the weakest material direction − π 4 and then merge in the middle of the specimen, according to the expectations. For ω = π 4 , this quite natural cracking pattern can no longer take place since the weakest material directions are now the vertical and the horizontal one. As a result, the two cracks propagate very closely to the vertical direction, and then merge along a direction close to horizontal. In both cases, no portion of the crack is aligned with the strongest material directions, hence no micro zig-zagging is observed. Once again the loss of symmetry in the final crack pattern for ω = π 4 suggests that alternative patterns are potentially possible upon slight perturbations [21].

Conclusions
Based on the family of Γ-convergent regularizations proposed by Focardi [16], we developed two variational phase-field models for fracture in anisotropic materials featuring two-fold and four-fold symmetric fracture toughness. Since both proposed models are of second order, as opposed to the previously available fourth-order models for four-fold symmetric fracture toughness, they do not require higher-continuity basis functions nor mixed variational principles for finite element discretization. However, the model for four-fold symmetry features a fracture energy functional that is non-convex with respect to the phase-field variable and thus requires a numerical solution algorithm able to cope robustly with non-convex minimization. In our numerical examples, we adopt the trust region method and assess the performance of the new models through several numerical examples simulating anisotropic fracture under anti-plane shear loading. The numerical results partly reflect the expectation that cracks propagate preferentially along the weakest material direction(s). In some cases, propagation is partially seen to occur macroscopically along strong directions, whereas microscopic zigzagging patterns (where 'microscopic' refers to a size scale of the order of the size of a finite element) are observed. Also, the results suggest an even higher sensitivity of the crack pattern to slight perturbations than already known for the isotropic fracture case [21].
Future work may include the following: • The proposed formulation for four-fold anisotropy should be numerically tested on more general loading conditions and possibly extended to the three-dimensional case. Also, its results should be thoroughly compared with those of the higher-order models for four-fold anisotropy available in the literature, such as the one in [33]; • The stochastic approach advocated in [21] to characterize non-uniqueness of the solution can be applied to the proposed anisotropic Foc-4 formulation, which seems to be even more sensitive to slight perturbations of the system than the formulations for isotropic fracture; • While this study concentrated on two-fold and four-fold anisotropy, other types of anisotropy, possibly featured in new artificial (e.g. 3D printed) materials can be in principle investigated with the same methodology.
A Existence and uniqueness of the phase-field solution for fixed displacement
Our findings are summarized in Table 4, whereas the derivations are outlined in the following. Table 4: Existence and uniqueness of the phase-field solution for fixed displacement for the isotropic Foc-4 and the anisotropic Foc-2 and Foc-4 formulations.
In the following, we first show that f satisfies a coercivity condition, and then prove that ξ → f (x, α, ξ) is convex for every (x, α) ∈ Ω × R, and (α, ξ) → f (x, α, ξ) is strictly convex for every x ∈ Ω. The first two properties are required for the existence, and the last one for the uniqueness of the solution.
The above results are summarized in Table 4. Notice that there is no correlation of this result with the convexity result for γ 2 in (6).
(63) conv_4a Here, the signs + and − are assumed to correspond to λ 1 and λ 2 , respectively. Owing to the estimates in (50), for every ω ∈ [0, π 2 ), the following holds As a result, λ 2 ≥ 0 iff τ ∈ [0, 1 3 ]. Using this along with (62), the existence of the solution to min α E can be concluded for the isotropic Foc-4 model (τ = 0) and for the anisotropic Foc-4 model when τ ∈ (0, 1 3 ]. Outside this range, the direct methods of the calculus of variations are insufficient to prove the existence result, and one may try to indirectly show existence by using, e.g., the Galerkin approximation [18]. In fact, this is what we do in Section 6.2 while performing numerical experiments for τ = 0.5, a value which lies outside the presumable solution existence range of [0, 1 3 ]. Finally, let us check the strict convexity of (α, ξ) → f (x, α, ξ) when τ ∈ [0, 1 3 ], required to establish the solution uniqueness in this range. The first eigenvalue of the corresponding Hessian reads λ 1 = g (α)Ψ(ε(u)) + 9G 0 b w α 2 , with g (α) = 56α 6 − 24α 2 , and the other two eigenvalues λ 2,3 coincide with those given by (63). The positivity of λ 1 in (x, α) ∈ Ω × R cannot be proven in general 4 . As a result, using the direct methods of the calculus of variations, we cannot conclude that min α E possesses a unique solution even for the solution existence range of τ ∈ [0, 1 3 ]. The above results are summarized in Table 4. Once again there is no correlation of this result with the convexity result for γ 4 in (7).