Spin-polarized transport in ferromagnetic multilayers: An unconditionally convergent FEM integrator.

We propose and analyze a decoupled time-marching scheme for the coupling of the Landau-Lifshitz-Gilbert equation with a quasilinear diffusion equation for the spin accumulation. This model describes the interplay of magnetization and electron spin accumulation in magnetic and nonmagnetic multilayer structures. Despite the strong nonlinearity of the overall PDE system, the proposed integrator requires only the solution of two linear systems per time-step. Unconditional convergence of the integrator towards weak solutions is proved.

Introduction

The interaction between electric current and magnetization in magnetic nanostructure devices and the control of this interaction have been realized through the prediction of the spin-transfer torque by Slonczewski and Berger  [1], [2]. The transfer of spin angular momentum between the spin-polarized electrical current and the local magnetization has been observed in various magnetic devices, such as metallic spin-valves systems, magnetic tunnel junctions, and magnetic domain walls in permalloy nanowires  [3], [4]. Based on these experiments, a number of technological applications have been proposed, e.g., STT-MRAMs, racetrack memories, and magnetic vortex oscillators  [5], [6].

The fundamental physics underlying these phenomena is understood as due to a spin torque that arises from the transfer of the spin angular momentum between conduction free electrons and magnetization. In the original works of Berger and Slonczewski  [1], [2], a homogeneous spin accumulation is assumed due to a current which flows through a first magnetic layer perpendicular to the interface into a second magnetic layer. The spin torque effect leads to an interaction between the spin-polarized current and the magnetization in the second layer. For magnetic multilayers it has been shown that a proper description of the magnetoresistance is essential to take into account the interplay between successive interfaces  [7], [8], [9]. In order to calculate the spin torque transfer, the spin transport properties have to be calculated far beyond the interface.

The original model of Berger and Slonczewski has been extended by taking into account the diffusion process of the spin accumulation by Shpiro et al. for one-dimensional systems  [10] and by García-Cervera and Wang  [11], [12] for three-dimensional systems. There, the overall system of PDEs (SDLLG) is a quasilinear diffusion equation for the evolution of the spin accumulation coupled to the Landau–Lifshitz–Gilbert equation (LLG) for the magnetization dynamics. Existence of global weak solutions to LLG goes back to  [13], while, in the same spirit, existence of global weak solutions to SDLLG is proved in  [11].

The reliable numerical integration of LLG (and, in particular, SDLLG) faces several challenges due to the nonuniqueness of weak solutions, the explicit nonlinearity, and an inherent nonconvex modulus constraint. Numerical approximation schemes for weak solutions of LLG are first proposed in  [14], [15]. First unconditional convergence results can be found in [16], [17], which consider the small-particle limit of LLG with exchange only. On the one hand, the integrator of  [16] relies on the midpoint rule and reduced integration, and thus has to solve one nonlinear system of equations per time-step. On the other hand, the tangent plane integrator of  [17], which extends the prior works  [14], [15], relies on a reformulation of LLG which is solved for the discrete time derivative. Each time-step consists of the solution of one linear system of equations plus nodal projection. It has been generalized to linear-implicit time integration and full effective field in  [18], [19].

Numerical integration of the coupling of LLG to other time-dependent PDEs has been analyzed in  [20], [21] for the full Maxwell equations (MLLG), in  [22], [23] for the eddy current formulation, and in  [24] for LLG with magnetostriction. While  [20] analyzes an extension of the midpoint scheme of  [16], the works  [21], [22], [23], [24] extend the tangent plane scheme from  [17], and emphasis is on the decoupling of the time-marching scheme in  [21], [23], [24].

In the models and works mentioned, e.g., MLLG, the coupling of LLG and Maxwell equations is weak in the sense that the magnetization of LLG only contributes to the right-hand side of the Maxwell system, while the magnetic field from the Maxwell equations gives a contribution to the effective field of LLG. In SDLLG the principal part of the differential operator of the spin diffusion equation depends nonlinearly on the magnetization. A first numerical integrator for SDLLG is proposed and empirically validated in  [12]. While this scheme appears to be unconditionally stable, the work does not prove convergence of the discrete solution towards a weak solution of SDLLG.

In our work, we extend the tangent plane integrator to SDLLG and prove unconditional convergence. Altogether, the contributions of the current work can be summarized as follows:

The proposed integrator is proven to converge (at least for a subsequence) towards a weak solution of SDLLG. This convergence is unconditional, i.e., there is no CFL-type coupling of the time and space discretizations. Despite the nonlinearity of SDLLG, each time-step requires only the solution of two successive linear systems, one for (the discrete time derivative of) the magnetization and one for the spin accumulation.

Our analysis thus provides, in particular, an alternate proof for the existence of (global) weak solutions of SDLLG, which has first been proved in  [11]. In addition to  [11], we prove that any weak limit of the proposed integrator satisfies an energy estimate similar to the theoretical behavior of (formal) strong solutions of SDLLG.

Unlike prior work on the tangent plane integrator, we adopt an idea from  [25] and show that the nodal projection step of the tangent plane scheme is not necessary. In particular and unlike the cited works, our analysis can therefore avoid a technical angle condition on the triangulations used. This result also transfers to the models and analysis of  [17], [18], [19], [21], [24], [23] and simplifies their (extended) tangent plane integrators.

Outline

The paper is organized as follows: In Section  2, we introduce and accurately describe the mathematical model, see (9) for the nondimensional formulation of SDLLG. In Section  3, we formulate a decoupled time-marching scheme (Algorithm 6) for the numerical integration of SDLLG and prove its well-posedness (Proposition 9). Section  4 contains the main result of our work (Theorem 12), which states unconditional convergence of the scheme towards weak solutions of SDLLG. Following  [11], weak solutions of SDLLG have finite energy. In Section  5, we prove that any weak limit obtained by the proposed numerical integrator shows the same energy behavior as formal strong solutions of SDLLG (Theorem 24). Numerical examples as well as the empirical validation of the proposed algorithm are postponed to a forthcoming paper  [26].

Notation

We use the standard notation  [27] for Lebesgue and Sobolev spaces and norms. For any domain , we denote the scalar product by for all . In the case of (spaces of) vector-valued functions, we use bold letters. For a sequence in a Banach space and , we write (resp.  ) in if the sequence converges strongly (resp. weakly) to in . Similarly, we write (resp. ) in if there exists a subsequence of which converges strongly (resp. weakly) to in . Throughout the paper, denotes a generic positive constant, independent of the discretization parameters, not necessarily the same at each occurrence. Alternatively, we write to abbreviate . Given , we denote by the tensor product defined by for all . By , we denote both the Frobenius norm of a matrix and the Euclidean norm of a vector. Since the meaning is clear from the argument, this does not lead to any ambiguity.

Model problem

In this section, we present the mathematical model, for which we introduce a nondimensional formulation, as well as the notion of a weak solution. We use physical units in the International System of Units (SI).

Physical background

We consider a magnetic multilayer. Let be polyhedral Lipschitz domains in , where corresponds to the volume occupied by the multilayer, and corresponds to the ferromagnetic part. A possible experimental setup is shown in Fig. 1. Given some finite time , we consider the time–space domains and .

Fig. 1

Schematic of a magnetic nanopillar structure (trilayer) consisting of two ferromagnetic films, and , separated by a nonmagnetic interlayer . The current is assumed to flow perpendicularly from to a bottom electrode connected to . In this case, and .

In micromagnetics, the quantity of interest is the magnetization , measured in ampere per meter (A/m). If the temperature is constant and far below from the Curie temperature of the ferromagnetic material, is a vector field of constant modulus , with being the saturation magnetization (in A/m). In the absence of spin currents, the dynamics of is described by the Landau–Lifshitz–Gilbert equation (LLG), which, in the so-called Gilbert form, reads Here, (radian per second per tesla) and (Newton per square ampere) are the gyromagnetic ratio and the permeability of vacuum, respectively, while is the nondimensional empiric Gilbert damping parameter. The effective field , measured in A/m, depends on and is proportional to the negative functional derivative of the total magnetic Gibbs free energy with respect to , i.e.,  In (2) the energy functional reads and consists of four terms, which correspond to the exchange energy, the anisotropy energy, Zeeman’s energy, and the magnetostatic energy, respectively. In (3), is the so-called exchange stiffness constant, measured in joule per meter (J/m), and is the anisotropic constant (in J/m3), while is a (nondimensional) smooth function, which takes into account the anisotropy of the ferromagnetic material. Moreover, is a given external field (in A/m), while refers to the magnetostatic potential, which is the unique solution of the full-space transmission problem Combining (2), (3), we obtain the following expression for the effective field: where denotes the stray field (in A/m).

The dynamics of the spin accumulation , measured in A/m, is described by the diffusion equation where is the diffusion coefficient (in m2/s), is the characteristic length of the spin-flip relaxation, and is related to the mean free path of an electron (both measured in m). The spin current , measured in A/s, is defined by where is the Bohr magneton, is the charge of the electron, and is the applied current density field (in A/m2), while the constants are the nondimensional spin polarization parameters of the magnetic layers. In (6) we denote by the matrix–vector product between the transpose of the Jacobian and , i.e., . In (5)–(6), it is implicitly assumed that in the nonmagnetic but conducting material .

To describe the dynamics of the magnetization, we take into account the interaction between the spin accumulation and the magnetization. Thus, we consider an augmented version of (1), namely where the constant in N/A2 is the strength of the interaction between the spin accumulation and the magnetization. Finally, to complete the setting, (5), (6), (7) are supplemented by initial conditions for some given initial states and with , and homogeneous Neumann boundary conditions

Nondimensional form of the problem

We introduce a nondimensional form of the system (5), (6), (7). We perform the substitution , with being the so-called (nondimensional) reduced time, and set . We rescale the spatial variable by , with being a characteristic length of the problem (measured in m), e.g., the intrinsic length scale . However, to simplify our notation, we write , and , instead of , and , respectively. We introduce the nondimensional vector unknowns , so that the modulus constraint becomes , and . Furthermore, we set and . With these substitutions, the nondimensional augmented form of LLG becomes where the effective field is given by with and , while the diffusion equation (5) reads To simplify our notation and without loss of generality, we assume that .

To sum up, we seek for with and such that Here, and are constants. For the diffusion coefficient , we assume that there exists a positive constant such that a.e. in . We also assume that and . Moreover, in (9a) we allow a more general effective field of the form where is a general time-independent field contribution. We emphasize that (10) in particular covers (8) with .

The constraint directly follows from the PDE formulation, provided in . Indeed, from (9a), we deduce that in .

Weak solution of the problem

Let be the dual space of and denote by the corresponding duality pairing, understood in the sense of the Gelfand triple . In view of the weak formulation of (9b), we consider the time-dependent bilinear form defined by for all and .

We recall from  [11, Definition 1] the notion of a weak solution of the SDLLG system (9), which extends the definition of weak solutions of LLG from  [13].

Let with a.e. in , and . The pair is called a weak solution of SDLLG if the following properties (i)–(iv) are satisfied:

with a.e. in and in the sense of traces,

and in the sense of traces,

for all , it holds

for almost all and all , it holds

If is a weak solution of SDLLG, then it holds and , cf., e.g.,  [27, Section 5.9.2, Theorems 2 and 3].

The following lemma highlights the parabolic nature of Eq. (9b).

The boundary term in (11b) is missing in  [11]. This error has recently been noticed and corrected, so that the overall result of  [11] remains valid  [28]. The present analysis provides an alternate proof for the existence of solutions and hence validity of the results of  [11], [28].

The time-dependent bilinear form is continuous and positive definite. Indeed, it holdsfor almost all .

The continuity directly follows from the regularity assumptions on the data, as a.e. in . As for the positive definiteness, we note As a consequence, since and , we get This establishes (12) and concludes the proof. □

Numerical algorithm

For the time discretization, we consider a uniform partition of the time interval with time-step size , i.e., for .

Given a sequence of functions , such that any is associated with the time-step , for we define the difference quotient . We consider the piecewise linear and the two piecewise constant time-approximations defined as follows: for and , we have Obviously, it holds for all .

For the spatial discretization, let be a shape-regular and (globally) quasi-uniform family of regular tetrahedral triangulations of , parameterized by the mesh size , where for all . By , we denote the restriction of to . We assume that is resolved, i.e.,  Let us denote by the standard finite element space of globally continuous and piecewise affine functions from to . Correspondingly, we also consider . By and , we denote the nodal interpolation operators onto these spaces. Since is resolved, these operators coincide on , i.e., for all . In particular, there is no ambiguity, if we denote both operators by . The set of nodes of the triangulation is denoted by .

We recall that, under the constraint , the strong form of (9a) can equivalently be stated as This formulation is used to construct the upcoming numerical scheme. Since (14) is linear in , the main idea is to introduce an additional free variable . To discretize , we introduce the discrete tangent space defined by for any . Moreover, we consider the set These sets reflect two main properties of and , namely the orthogonality and the unit-length constraint .

Let . We consider the nodal projection map defined by for all and . A simple argument based on the elementwise use of barycentric coordinates shows that for all . Moreover, we have the estimate where the constant depends only on the shape-regularity of the triangulation, cf., e.g.,  [25, Lemma 2.2]. With an additional angle condition on , it is well known that (15) holds even with , cf.  [29].

Let and be suitable approximations of the initial conditions. Moreover, we consider a numerical realization of , which is assumed to fulfill a certain set of properties, see (H2)–(H3) below. This allows us to include the approximation errors, e.g., those which arise from the numerical computation of the stray field, into the overall convergence analysis. For ease of presentation, we assume that and are continuous in time, i.e., and , so that the expressions and are meaningful for all . It is even possible to replace and by some numerical approximation and as long as some weak convergence properties are fulfilled, cf.  [19].

Analogously to what we have done in Section  2 for the continuous problem, for we define the bilinear form by for all . For the numerical integration of the SDLLG system (9), we propose the following algorithm. The overall system (9) is a nonlinearly coupled system of a linear diffusion equation for with the nonlinear LLG equation for . However, our scheme only requires the solution of two linear systems per time-step, since the treatment of the micromagnetic part and the spin diffusion part is completely decoupled for the time-integration. This greatly simplifies an actual numerical implementation as well as the possible preconditioning of iterative solvers. The following result follows from standard scaling arguments. The following proposition states that the above algorithm is well defined, cf.  [25, Propositions 3.1 and 4.1] for corresponding results in the frame of harmonic maps and the harmonic map heat flow.

Input: , parameter .

For all iterate: Output: Sequence of discrete functions .

compute such that for all ;

define by

compute such that for all .

Unlike this work, earlier results on the tangent plane integrator  [17], [18], [19], [20], [21], [22], [23], [24] define in (16b). Unconditional convergence in the sense of Theorem 12 can then be achieved with an additional angle condition on the triangulation , which ensures (15) with . This assumption is avoided in the present work.

Let be a quasi-uniform family of triangulations of and . Then,The constant depends only on , but is independent of the mesh size . □

Algorithm  6   is well defined in the following sense: for each time-step , there exists a unique solution . Moreover, it holdsas well aswhere the constant depends only on the shape-regularity of , but is independent of and .

Let . For step (i) of the algorithm, it is straightforward to show that problem (16a) is characterized by a positive definite bilinear form. Unique solvability thus follows from linearity and finite space dimension. Step (ii) is clearly well defined. For all , the nodewise orthogonality from proves Since , mathematical induction proves This proves (17). The norm equivalence from Lemma 8 in the case yields This establishes (18). For step (iii), we use the same argument as for step (i). Due to (17), the nodewise projections in (16c) are well defined. Let be the bilinear form associated to problem (16c), i.e.,  Since , we see It follows that As and , is positive definite and problem (16c) is thus well posed. □

Convergence analysis

In this section, we consider the convergence properties of Algorithm 6 and show that it is indeed unconditionally convergent towards a weak solution of SDLLG in the sense of Definition 2. We emphasize that the proof is constructive in the sense that it even shows existence of weak solutions. We start by collecting some general assumptions:

The discrete initial data and satisfy

The general field contribution is bounded, i.e.,  with a constant which depends only on .

It holds for any sequence in .

Usual stray field discretizations by hybrid FEM–BEM methods, e.g., the Fredkin–Koehler approach from  [30], or FEM–BEM coupling methods satisfy (H2)–(H3), see  [19].

From now on, we consider the time-approximations defined by (13). The next theorem is the main result of this work.

For a discrete operator , assumption (H2) can be relaxed to Within this setting, and with an appropriate modification of assumption (H3), the hybrid FEM–BEM method from  [31] for the computation of the stray field can also be included into our analysis. Then, the proof of Proposition 19 becomes more technical, but the assertion remains true. We refer to the argument of  [19] which can be adapted accordingly.

Let be a shape-regular and quasi-uniform family of triangulations.

Suppose and that assumptions  (H1)–(H2)  are satisfied.

Then, there exist and such that

In addition to the above, let assumption  (H3)  be satisfied. Then, it holdswhere is a weak solution of SDLLG.

The proof of Theorem 12 will roughly be done in three steps, namely For the sake of readability, we split our argument into several lemmata.

In particular, Theorem 12 yields existence of weak solutions, and each accumulation point of is a weak solution of SDLLG in the sense of Definition 2.

boundedness of the discrete quantities and energies,

existence of weakly convergent subsequences via compactness,

identification of the limits with weak solutions of SDLLG.

To start with, we recall the following result, which states a well-known and simple algebraic trick which often simplifies the computation and the estimation of sums. The first ingredient for step (i) is the following proposition.

Abel’s Summation by Parts

Let be a vector space endowed with a symmetric bilinear form . Given an integer , let . Then, it holds □

Let and suppose that the assumptions of   Theorem  12(a)  are satisfied. Then, the discrete functions obtained through   Algorithm  6   fulfillThe constant depends only on the data, but is in particular independent of the discretization parameters and .

Let . For (16c), we choose as test function. After multiplication by , we obtain Since and , it follows that cf. the proof of Lemma 5. Summing up over , and exploiting Abel’s summation by parts from Lemma 14 for the term , we get Exploiting on the left-hand side, the Cauchy–Schwarz inequality and the Young inequality on the right-hand side, we obtain, for any choice of , Here the constant is the stability constant of the trace operator. It follows that If we choose , then all the coefficients on the left-hand side are positive. From (H1) and the regularity of , we know that the right-hand side is uniformly bounded with respect to and . This yields the estimate (19). □

Under the assumptions of   Proposition  15, the sequences and are uniformly bounded in and in , i.e.,where the constant depends only on the data, but is in particular independent of the discretization parameters and .

Let be the -orthogonal projection onto , i.e.,  Since is quasi-uniform, it is well known that is stable in , i.e.,  We also refer to  [32], [33] for recent results on -stability on locally refined meshes. With this, we obtain uniform boundedness of .

The result follows from the boundedness of the discrete functions from Proposition 15. □

The sequence is uniformly bounded in , i.e.,where the constant depends only on the data, but is in particular independent of the discretization parameters and .

We derive the corresponding estimates for the discrete quantities .

Let , , and . From (16c) and the -stability (20) of , we get Dividing by and taking the supremum over , we obtain Squaring, integrating over , and summing up over , we get The boundedness from Proposition 15 thus yields (21). □

Let . The discrete functions obtained through   Algorithm  6   fulfill

We test (16a) with to get Exploiting the vector identity with the choice and , and taking into account (16b), we obtain (22). □

Suppose that the assumptions of   Theorem  12(a)  are satisfied. Then, there exists such that for all time-step sizes and the discrete functions obtained through   Algorithm  6   fulfillThe constant depends only on the data and , but is otherwise independent of the discretization parameters and .

Let . From Lemma 18, multiplying (22) by , summing up over and exploiting the telescopic sum, we obtain The Cauchy–Schwarz inequality and the Young inequality, together with assumption (H2), yield for any From Proposition 9, we deduce where the constant depends only on and . We thus obtain Note that . If we choose , for all the coefficients on the left-hand side are positive. From the regularity of , assumption (H1), and the boundedness from Proposition 15, we know that the right-hand side is uniformly bounded. This yields the estimate (23). □

Under the assumptions of   Proposition  19, and if , the sequences and are uniformly bounded. In particular, it holdswhere the constant depends only on the data and , but is independent of the discretization parameters and .

We can now proceed with step (ii) of the proof and conclude the existence of weakly convergent subsequences.

The result follows from the boundedness of the discrete functions from Proposition 9, Proposition 19, and from (15). □

Suppose that the assumptions of   Theorem  12(a)  are satisfied. Then, there exist and , with a.e. in , such that there holdsfor . Moreover, there exists one subsequence for which   (24)   holds simultaneously.

The boundedness results from Corollary 20, in combination with the Eberlein–Šmulian theorem, allow us to extract weakly convergent subsequences of and . Let be such that in . From the continuous inclusions and the compact embedding , we deduce With , we can identify the limits of the subsequences of and . As , it clearly holds that a.e. in .

We now prove that the limiting function satisfies the unit-length constraint. First, we observe that for . For and , a standard interpolation estimate for the piecewise linear function yields From Proposition 19, we obtain For all , Proposition 9 and the discrete norm equivalence of Lemma 8 with yield Then, from Proposition 19, we deduce Combining (26)–(27), the triangle inequality thus yields that in for , whence a.e. in follows from (25).

For with , it holds that Due to (17), this yields for all whence by virtue of Proposition 19 This implies in as . Since , we have even in as well as in .

From Corollary 16, we similarly deduce the existence of weakly convergent subsequences of and . Due to Proposition 15, the quantity is bounded. This allows to identify the weak limits, since Finally, from Proposition 17, we deduce the existence of a weakly convergent subsequence of , and it is easy to see that its limit is precisely , cf.  [27, Section 7.1.2, Theorem 3]. This establishes (24e)–(24f) and thus concludes the proof. □

We have collected all the ingredients for the proof of our main theorem.

As the constants which guarantee the boundedness of Proposition 15, Corollary 16 and Proposition 17, are independent of , we deduce that .

The result of part (a) follows directly from Proposition 21. To conclude the proof of part (b), it remains to identify the limiting functions with a weak solution of SDLLG in the sense of Definition 2.

To check (11a), we essentially proceed as in  [17]. Let . For , we test (16a) with respect to , with being the nodal interpolation operator onto . Multiplication with and summation over yield where for all . From the well-known approximation properties of and the boundedness of from Proposition 19 for , we deduce Passing to the limit for , we obtain In the latter, we have used the convergence properties from Proposition 21, assumption (H3) for the general field contribution, as well as and in .

Direct calculations and standard properties of the cross product yield the identities from which, by density, we deduce (11a).

To check (11b), let . Given , let . In (16c) we choose the test function . Integration in time over and summation over yield where for all . Passing (28) to the limit for , due to the convergence properties stated in Proposition 21, in combination with the standard approximation properties of , we deduce By density, this is also true for all . Hence in particular for each and a.e.  we have (11b).

Since in with , and in with , assumption (H1) allows to deduce and in the sense of traces. □

Energy estimate

In this section, we exploit our constructive convergence proof to derive an energy estimate for weak solutions of SDLLG, which is also meaningful from a physical point of view. The total magnetic Gibbs free energy from (3) is strongly related to the standard form (1) of LLG and does not take into account the interaction between the magnetization and spin accumulation. As we are dealing with the augmented form (7) of LLG, we extend (3) and define the free energy of the system by This definition is in agreement with (7), since in this case it is easy to see that where the effective field is given by (4). A simple formal computation shows that strong solutions to (7) satisfy Neglecting the spin accumulation term and assuming that is constant in time, Eq. (30) reduces to which reveals the well-known dissipative behavior of solutions to the standard form (1) of LLG.

The main aim of this section is to prove a property corresponding to (30) in the context of weak solutions. To this end, we move to the nondimensional framework introduced in Section  2 and consider the following assumptions:

the operator from (10) is linear, self-adjoint, and bounded;

it holds for any sequence in , which is slightly stronger than (H3);

the applied field belongs to .

Up to an additive constant, the nondimensional counterpart of (29) reads The following theorem proves an energy estimate which generalizes (30) to weak solutions.

For some fixed easy axis and the corresponding uniaxial anisotropy density function , the operator satisfies (H4). Moreover, all the stray field discretizations mentioned in  Remark 10 satisfy (H5), see  [19]. The operator from (31) is even well defined and bounded as operator for all , see  [34]. Unlike  [18], the proof of our energy estimate, see Theorem 24, avoids this additional regularity, but only relies on the energy setting .

Suppose that assumptions  (H1)–(H2)  and  (H4)–(H6)  are satisfied. Let be a weak solution of SDLLG obtained as a weak limit of the finite element solutions from   Algorithm  6   for . Then, the energy functional from   (32)   satisfiesfor almost all .

Given , let such that . Let . From Lemma 18, we get By definition (16b), it holds . We thus obtain Analogously, we see that . Since is linear and self-adjoint, (16b) also reveals Altogether, we thus obtain Since is a bounded operator, it follows that Summation over and the boundedness from Proposition 19 yield The available convergence results on , and , as well as assumption (H5), allow us to employ standard arguments with weakly lower semicontinuity for the limit . This concludes the proof of (33). □