{'fraction_non_alphanumeric': 0.05167450677525802, 'fraction_numerical': 0.018324791125465142, 'mean_word_length': 4.344840525328331, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Revisit of The Energy Quadratization Method with A Relaxation Technique

This letter revisits the energy quadratization (EQ) method by introducing a novel and essential relaxation technique to improve its accuracy and stability. The EQ method has witnessed significant popularity in the past few years. Though acknowledging its effectiveness in designing energy-stable schemes for thermodynamically consistent models, the primary known drawback is apparent, i.e., its preserves a"modified"energy law represented by auxiliary variables instead of the original variables. Truncation errors are introduced during numerical calculations so that the numerical solutions of the auxiliary variables are no longer equivalent to their original continuous definitions. Even though the"modified"energy dissipation law is preserved, the original energy dissipation law is not guaranteed. In this paper, we overcome this issue by introducing a relaxation technique. The computational cost of this extra technique is negligible compared with the baseline EQ method. Meanwhile, the relaxed-EQ method holds all the baseline EQ method's good properties, such as linearity and unconditionally energy stability. Then we apply the relaxed-EQ method to several widely-used phase field models to highlight its effectiveness.


Background
In this letter, we focus on the PDE models that respect the thermodynamical laws. They are usually named thermodynamically consistent models. In particular, when the temperature deviation is considered, the PDE model's entropy shall be non-decreasing in time. When the temperature deviation is ignored, the PDE model's Helmholtz free energy shall be non-increasing in time. These thermodynamically consistent models are broadly used to investigate problems in various fields, taking material science, life science, and mechanical engineering as examples.
To illustrate the idea of the EQ method, we start with a simple example, the Allen-Cahn (AC) equation with periodic boundary condition, given as ∂ t φ(x,t) = ε 2 ∆φ(x,t)−φ(x,t) 3 +φ(x,t), (x,t) ∈ Ω×(0,T]. (1.1) Define the free energy To design numerical schemes that preserve the energy dissipation law for the AC equation, the EQ method introduces a novel perspective, namely, instead of solving the original AC equation in (1.1), it solves an equivalent problem. By introducing the auxiliary variable q(x,t) := h(φ) = √ 2 2 (φ 2 −1), and the notation g(φ) := ∂q ∂φ = √ 2φ, the equivalent problem is formulated as where the modified free energy is defined as E = Ω , as well as their energy dissipation laws.
However, the numerical solutions for q and the numerical results of h(φ) are not necessarily equal any more (in fact, they are not) after temporal discretization. Hence a discrete energy dissipation law in (1.4) is no longer equivalent with a discrete energy dissipation law in (1.2), though we shall recognize that they are consistent with the numerical schemes' order of temporal accuracy.
With this in mind, the primary goal of this letter is to overcome this inconsistency issue. And our central idea is to introduce a relaxation technique for the numerical solutions of q, to penalize the numerical deviation from the original definition of q := h(φ).

A generic model formulation
Consider a domain Ω, and denote the thermodynamic variable as φ. We introduce the general formulation for a thermodynamically consistent reversible-irreversible PDE models as where B is a trace operator, and G := G a +G s is the mobility operator, G s is symmetric and positive semi-definite to ensure thermodynamically consistency, G a is skew-symmetric, and δE δφ is the variational derivative of E , known as the chemical potential. Then, the triplet (φ,G,E ) uniquely defines a thermodynamically consistent model. One intrinsic property of (2.1) owing to the thermodynamical consistency is the energy dissipation law where the inner product is defined by ( f ,g) = Ω f gdx,∀ f ,g ∈ L 2 (Ω), andĖ sur f is due to the boundary contribution, and g b is the boundary integrand. When G a = 0, (2.1) is a purely dissipative system; while G s = 0, it is a purely dispersive system.Ė sur f vanishes only for suitable boundary conditions, which include periodic and certain physical boundary conditions.

Model reformulation with the energy quadratization method
Now, we illustrate the energy quadratization (EQ) idea on solving (2.1). Denote the total energy as E (φ) = Ω edx, with e the energy density function. We denote L 0 as a linear operator that can be separated from e. Introduce the auxiliary variable where A 0 > 0 such that q is a well defined real variable. Then we rewrite the energy as With the EQ approach above, we transform the free energy density into a quadratic one by introducing an auxiliary variable to "remove" the quadratic gradient term from the energy density. Assuming q = q(φ,∇φ) and denoting g(φ) = ∂q ∂φ and G(∇φ) = ∂q ∂∇φ , we reformulate (2.1) into an equivalent form with the consistent initial condition q| t=0 = 2 e− 1 2 |L . Now, instead of dealing with (2.1) directly, we develop structure-preserving schemes for (2.5).
The advantage of using model (2.5) over model (2.1) is that the energy density is transformed into a quadratic one in (2.5).
where L is a linear operator, and N (Ψ) is the mobility operator as

Model examples
To demonstrate this methodology's generality, we introduce a few widely-used PDE models that can be cast as special cases.
Allen Cahn equation. First of all, let us consider the Allen-Cahn equation It can be formulated in an energy variation form of (2.1) with G = 1 and F = Ω Then the model is reformulated into the quadratization form of (2.5) as with the consistent initial condition q| t=0 = h(φ| t=0 ). Its matrix formulation of (2.6) can also be easily obtained, by noticing L 0 = −ε 2 ∆+γ 0 . Cahn-Hilliard (CH) equation. In the second example, we consider the Cahn-Hilliard equation given as where f (φ) is the bulk potential. In this paper, we choose the double-well potential f (φ) = 1 4 (φ 2 − 1) 2 . Other potentials, such as the Flory-Huggins potential can also be used. It can be reformulated in an energy variational form of (2.1) with the notations G = −∆, and F = Ω 2φ, leading us to the quadratized form of (2.5) as with the consistent initial condition q| t=0 = h(φ| t=0 ). Molecular beam epitaxy (MBE) model. Then, we focus on the molecular beam epitaxy (MBE) model with slope section. The model reads Similarly, this model can be written in an energy variation form of (2.1) as G=M and F= Ω 2∇φ, the MBE model can be written in the quadratized form as with a consistent initial condition for q. Phase field crystal (PFC) equation. Next, we investigate the phase field crystal (PFC) equation which reads as with a 0 and b 0 are model parameters. The PFC model can also be written in an energy variational form of (2.1) with G = ∆ , and F = Ω , we can rewrite the PFC model into the quadratizated form of (2.5) as 3 The relaxed energy quadratization method

Baseline EQ schemes
We recall the classical EQ numerical algorithms to solve the generalized thermodynamically consistent problems. Mainly, we discretize all the linear terms implicitly and all the nonlinear terms explicitly. We will end up with some energy stable and linear numerical algorithms. 2 )LΨ n+ 1 2 , with the notation Ψ n+ 1 2 = 3 2 Ψ n − 1 2 Ψ n−1 . Scheme 3.2 (BDF2 Scheme). After we calculate Ψ n−1 and Ψ n , we can update Ψ n+1 via the following scheme It could be easily verified that the EQ schemes are uniquely solvable and unconditionally energy stable [10]. Meanwhile, the known issue is that energy is modified, since q does not necessarily equal to its original definition of (2.3) after discretization, i.e. q n+1 = h(φ n+1 ) in general.

Relaxed EQ schemes
To overcome this inconsistency issue of the modified energy and the original energy in the baseline EQ method above, we come up with the novel relaxed EQ schemes as follows, as an analogy to the two baseline EQ schemes. • Step 2, update the baseline solutions with a relaxation step. We update q n+1 as where ξ 0 is calculated from with the coefficient a, b and c given by Here η ∈ [0,1] is an artificial parameter that can be assigned. Remark 3.1. Here we explain the idea for the relaxation step. ξ 0 is a solution for the optimization problem: It is not hard to show that ξ = 1 is in the feasible domain, since a+b+c ≤ 0. In addition with a > 0, we can easily have ∆ = b 2 −4ac ≥ 0. Then we have two real roots for the quadratic equation aξ 2 +bξ +c =0 that are given by ξ = −b± ≤ 1, so that we obtain the solution for (3.2) as Remark 3.2. We emphasize that η ∈ [0,1] in (3.2) controls the relaxation strength. When η = 0, a minimum relaxation is introduced. When η = 1 a maximum relaxation is introduced. Remark 3.3. The Scheme 3.3 is second-order accurate in time. Notice the fact,q n+1 =q(t n+1 )+O(δt 2 ) and h(φ n+1 ) = q(t n+1 )+O(δt 2 ). Thus, ∀ξ ∈ [0,1], ξq n+1 +(1−ξ)h(φ n+1 ) = q(t n+1 )+O(δt 2 ). I.e., the relaxation step doesn't affect the order of accuracy for the proposed schemes.
Remark 3.4. We point out the following connections. When ξ 0 = 1, the Scheme 3.3 reduces to the baseline CN scheme 3.1. When ξ 0 ∈(0,1), the Scheme 3.3 relaxes the solution q n+1 to be closer to its original definition of (2.3) in the continuous level. And when ξ 0 =0, the solution q n+1 is consistent to its original definition of (2.3). In practice, the artificial parameter η ∈[0,1] is used to control how much relaxation is added.
Theorem 3.1. The Scheme 3.3 is unconditionally energy stable.
Proof. As a matter of fact, if we take inner product of (3.1) with δtLΨ n+ 1 2 , we will have Also, from (3.2), we know Adding two equations (3.5) and (3.6) together gives us the energy inequality Similarly, we can propose the relaxed BDF2 scheme as follows.
Scheme 3.4 (Relaxed BDF2 Scheme). After we calculate Ψ n−1 and Ψ n , we can update Ψ n+1 via the following two step schemes: • Step1, we obtain the intermediate solutions via the baseline EQ method. Specifically, we • Step2, we update the baseline solutions via a relaxation step. We update q n+1 by Here η ∈ [0,1] is an artificial parameter that can be manually assigned.
Remark 3.5. In a similar manner, the readers shall notice the extra relaxation step (3.9) is cheapand-easy to compute. We comment on the idea for the relaxation step. ξ is a solution of the optimization problem where η ∈ [0,1] is an artificial parameter that can be manually assigned.
Remark 3.6. In practice, it can be simplified as solving the following algebra optimization problem 1] ξ, subject to aξ 2 +bξ +c ≤ 0, (3.12) where the coefficients for the quadratic equation are By a similar argument, we have Theorem 3.2. The Scheme 3.4 is unconditionally energy stable.
Proof. The proof is similar with the proof for Scheme 3.3. In fact, if we take inner product of (3.8) with δtΨ n+1 , we will have Meanwhile, we know from (3.9) that Adding the equations (3.14) and (3.15) together, we will have This completes the proof.

Numerical Results
Next, we test this novel relaxed-EQ method on several specific thermodynamically consistent models. We consider problems with periodic boundary conditions for simplicity, though we emphasize our approach has no specific restriction on the boundary conditions. For spatial discretization, we use the Fourier Pseudo-spectral method. Due to space limitation, we only test the CN schemes, i.e. Scheme 3.1 and Scheme 3.3 and choose η = 1.
For the AC equation and the CH equation, we consider a square domain Ω = [0,1] 2 with seven disks in the 2D domain as an initial condition, similar with the problem in [1]. It is known that the disks will shrink and eventually disappear for the dynamics driven by the AC equation since the AC equation doesn't preserve the total volume. Meanwhile, the Ostwald ripening effect will occur for the CH equation dynamics since the CH equation maintains the total volume while minimizing the interfaces. Numerical results for the two dynamics are summarized in Figure 4.1(a) and 4.1(b), respectively, where we observe the expected dynamics.  One particular focus of this letter is how the relaxed-EQ method performs compared with the baseline EQ method. We solve the problem in Figure 4.1 using both the EQ and relaxed-EQ methods with 128 2 meshes and a large time step, namely, δt = 0.75 for the AC equation, and δt = 0.005 for the CH equation. The numerical energies using both approaches are summarized in Figure 4.2(a) and 4.2(b) respectively for the AC and CH equations. Notice both algorithms perform well when the time step is small enough. When the time step is large, the relaxed-EQ method outperforms the baseline EQ method.
Furthermore, we also use both methods to solve the MBE equation. We use classical benchmark problems from [5,8], with the same initial conditions and parameters. Due to space limitations, the numerical results are not shown. Instead, the energy evolutions using both numerical methods with a time step δt = 0.001 are summarized in Figure 4.2(c). We observe that the relaxed-EQ method provides more accurate results than the EQ method for solving the MBE model as well.
In the last example, we solve the PFC model using the relaxed-EQ scheme to further highlight  its power for simulating long-time dynamics. We set up the initial conditions as

Conclusion
This letter reports a relaxation technique to overcome the inconsistency issue of the modified free energy and the original energy in the energy quadratization (EQ) method. Through theorems and numerical examples, it is noticeable that the relaxed-EQ method always outperforms the baseline EQ method. In particular, it overcomes the EQ method's primary issue that the auxiliary variable deviates from its original continuous definition after discretization as numerical errors are introduced and accumulated. In the relaxed-EQ method, the numerical solution of the auxiliary variable is relaxed towards its continuous definition in each time step with negligible computational cost. Overall, this letter's main message is that the relaxation step shall be applied to all available EQ schemes in literature with its easy-to-use and accuracy-improving nature.