Investigation of Viscoelastic-Creep and Mechanical-Hysteresis Behaviors of Hydrostatically Stressed Crystal Using the Phase Field Crystal Method

The phase field crystal (PFC) method is a density-functional-type model with atomistic resolution and operating on diffusive time scales which has been proven to be an efficient tool for predicting numerous material phenomena. In this work, we first propose a method to predict viscoelastic-creep and mechanical-hysteresis behaviors in a body-centered-cubic (BCC) solid using a PFC method that is incorporated with a pressure-controlled dynamic equation which enables convenient control of deformation by specifying external pressure. To achieve our objective, we use constant pressure for the viscoelastic-creep study and sinusoidal pressure oscillation for the mechanical-hysteresis study. The parametric studies show that the relaxation time in the viscoelasticcreep phenomena is proportional to temperature. Also, mechanical-hysteresis behavior and the complex moduli predicted by the model are consistent with those of the standard linear solid model in a low-frequency pressure oscillation. Moreover, the impact of temperature on complex moduli is also investigated within the solid-stabilizing range. These results qualitatively agree with experimental and theoretical observations reported in the previous literature. We believe that our work should contribute to extending the capability of the PFCmethod to investigate the deformation problemwhen the externally applied pressure is required.


Introduction
The phase field crystal (PFC) method has emerged as a computational model with atomistic resolution and diffusive time scale. The method has an advantage over molecular dynamics (MD) in terms of the time scale that is not restricted by lattice vibration time; this is due to the specification of the order parameter as the local-time-average atomic number density and the evolution of the order parameter through dissipative dynamics. The method also has advantage over the phase field (PF) method in terms of atomic resolution and self-consistency; this benefit stems from the free energy functional that can be minimized by the periodic order-parameter field whereas the free energy functional of the PF method is formulated by a spatially uniform field which diminishes numerous physical features that occurs due to the periodicity of crystalline phases, e.g., spatial symmetry, elastic and plastic interaction, nucleation, multiple crystal orientations, and motion of dislocations [1][2][3]. The PFC method was first introduced by Elder et al. [1,2] and has been extensively used to study material phenomena and behavior ranging from phase transformation [4,5] and topological defect dynamics (vacancy [6], grain boundary [7][8][9], dislocation [3,10,11], and crack [12]) to elasticity and elastic constants [13][14][15] and plasticity [11,16]. Moreover, the PFC method was also interpreted and derived according to the classical density functional theory (CDFT) of the freezing point of view [4]. This derivation provided an additional field variable which extended the capability of the PFC method to investigate the material phenomena in a more complex situation such as phase transformation [4,17] and segregation-induced grain-boundary premelting [18] in binary alloys.
Although there were numerous successful PFC studies on the mechanical behavior of materials, such as elasticity, plasticity, and phase transformation in pure materials and binary alloys, the investigation of viscoelastic behavior and the related properties predicted by the PFC model is still limited and not thoroughly explored. Thus, it is of interest to study the capability of the PFC model to predict these phenomena. Viscoelasticity is the behavior with both viscous and elastic characteristics; thus, materials with viscoelastic behavior can return to their original configuration when unloaded but do so in a time-dependent manner. The viscous property also leads to the fact that material response depends on the rate at which it is deformed [19,20]. Viscoelastic effect has played significant roles in many materials ranging from amorphous polymers [21], semicrystalline polymers [22], biopolymers [23], and bitumen materials [24] to metals at very high temperature [25]. Furthermore, viscoelasticity governs many engineering applications in various fields including damping material design for noise reduction [26,27] and vibration control in engineering structure [28], shapememory ceramics and polymers [29], piezoelectric materials [30], self-healing materials [31], viscoelastic gels in surgery application [32], low-loss materials [33], and nanoscale resonators [34,35].
In quantifying the viscoelastic effect, there are three types of studies: viscoelastic creep, stress relaxation, and mechanical hysteresis. In this work, we not only focus on the viscoelastic creep but also extend our study to mechanicalhysteresis behavior and its related parameters in a bodycentered-cubic (BCC) solid (3D case) predicted by the PFC method, unlike our previous works [36,37] which were solely limited to viscoelastic creep in a one-dimensional (1D case) crystal study; moreover, the mechanical-hysteresis behavior and its related parameters were not addressed in those works. Since the studies of viscoelastic creep and mechanical hysteresis involve measuring the strain response from an applied stress, we also need a method to control the applied stress, and thus, we employ the PFC model that is incorporated with our previously proposed pressure controlled dynamic (PCD) equation [37]. The use of the PCD equation enables convenient control of deformation using external pressure which allows studies of experimentally observed phenomena involving a system under different types of applied pressure; we will refer to this particular PFC method as the PFC-PCD model. The PCD equation is similar to the dynamic equation introduced by Kocher and Provatas [38]. They used their PCD equation to control internal pressure by specifying external pressure to investigate the full spectrum of solidliquid-vapor transitions; unfortunately, the viscoelasticcreep and mechanical-hysteresis behaviors which can be produced from that equation were not the focus area and, of course, not addressed in their study [38]. The difference between two PCD versions is due to the different derivation process in which our version is established regarding the thermodynamics of a hydrostatically stressed crystal solid proposed by Larché and Cahn [39] and classical irreversible thermodynamic frameworks [40][41][42]. Moreover, the deviation of viscoelastic-creep behavior in 1D crystal between two PCD versions was also investigated in our previous study [37]. The advantage of the PCD equation is that it enables the external pressure to be specified as an input to the simulation, which allows the system volume, or the strain response, to be measured; this is in contrast to the conventional PFC simulations where volume, or grid spacing, is an independent variable [43].
For the investigation of viscoelastic behavior, we simulate a hydrostatic stress via a BCC solid subjected to two types of applied pressure: constant pressure for the viscoelastic-creep study and sinusoidal pressure oscillation for the mechanical-hysteresis study. We find that, in both cases, the solid exhibits delayed strain responses from the applied pressure which is indicative of viscoelasticity, and the degree of viscoelasticity can be controlled by the parameter in the PCD equation. Furthermore, the strain responses from the simulations can be compared with functional forms of the solutions to the standard linear solid model to extract viscoelastic quantities and properties. These properties include relaxation time of viscoelastic creep and the complex moduli from the hysteresis. The parametric studies show that the relaxation time is proportional to temperature while the magnitude of deformation at steady-state condition is proportional to both temperature and atomic density. Also, the dynamic behavior, mechanical hysteresis, and complex moduli predicted by the PFC-PCD equation under sinusoidal pressure oscillation are consistent with those by a standard linear solid model at a certain frequency range. However, at high frequency range, the behavior exhibited by the PFC-PCD model differs from that by the standard linear solid model. The reason behind this difference is thoroughly analyzed by conducting Taylor's expansion on PCD equation in Appendix C. Moreover, the impact of temperature on complex moduli is also investigated, and the results well agree with experimental and theoretical observation. We are positive that our work should provide the extension of capability for the PFC method to investigate the deformation problem in a real crystalline solid when externally applied pressure is necessary. This paper is organized as follows. In Section II, we provide the background information on the PFC method, the pressure-control dynamic equation, and the standard linear solid model. In Section III, the method of simulation is presented. The results and discussions are provided in Section IV, and we conclude the work with the summary in Section V.

Background
2.1. PFC Method. The PFC method is characterized by its free energy functional and the dynamic equation. The simplest free energy functional was introduced by Elder and coworkers [1,2]: Advances in Mathematical Physics wðρÞ is the free energy density; a t , g t , λ, and q 0 are fitting parameters; and ρ is the atomic number density field. This PFC free energy functional can be minimized by a constant ρ profile, representing a liquid phase, and a periodic ρ profile, representing a crystalline solid phase. The density periodic field ϕ can be expressed in terms of Fourier expansion of the form where G j is a reciprocal lattice vector (RLV), A j is an amplitude of the density wave corresponding to G j , ρ is an average atomic number density, and c.c. denotes a complex conjugate. The minimization of F essentially results in the truncation of the summation in Equation (2) to terms with small values of jG j j as the terms with higher values of jG j j have vanishingly small A j . Consequently, this particular PFC model (Equation (1)) favors crystal structures that can be constructed from a single set of RLVs such as stripe, hexagonal, and BCC structures in one, two, and three dimensions, respectively. To this end, one can construct the so-called "one-mode" approximation of the BCC crystal by limiting the terms in Equation (2) to those corresponding to h110i RLVs or jG j j = 2π√2/L a , where L a is the lattice parameter. By assuming that A j of the terms with the same jG j j are equivalent, one obtains ρ one = ρ + 4A s cos q a r 1 ð Þcos q a r 2 ð Þ+ cos q a r 1 ð Þcos q a r 3 ð Þ ½ + cos q a r 2 ð Þcos q a r 3 ð Þ, where q a = 2π/L a and A s is the density-wave amplitude.
To simplify the expressions in Equation (1), one can nondimensionalize the variables by using the following substitutions [2]:r where the quantities with tildes are nondimensional and d = 3 is the dimensionality of the problem. The expressions in Equation (1) simplify tõ where ϵ can be interpreted as the degree of undercooling which is inversely proportional to temperature. The onemode approximation of the BCC crystal becomes whereq a = q a /q 0 = 2π/ðL a q 0 Þ = 2π/L a . The value ofq a can be determined from minimization of the free energy density [44] leading toq a = 1/ ffiffi ffi 2 p and, consequently,L a = 2π ffiffi ffi 2 p .
One can also approximate the ranges of∈ and ρ where the BCC solid is stable. This is achieved by calculating the free energy density of the solids using the one-mode approximations for stripe, hexagonal, and BCC structures and the free energy of the liquid using a uniform density. Then, the phase diagram can be determined from the common tangent construction, and the ranges of∈ and ρ where the BCC solid is stable can be obtained [45].
The evolution equation forρ is the Cahn-Hilliard type which involves dissipative dynamics and mass conservation [2].
whereL μ is the mobility coefficient andμ = δF /δρ is the diffusion potential. This equation can be used to simulate evolution ofρ under specified (or fixed) temperature (from the value of∈), mass (from the value of ρ), and volume (from value of grid spacing); typical simulations under this condition are solidification and microstructural evolution. One can also vary the grid spacing to simulate the evolution of the density profile under specified timedependent deformation [43]. In the context of the viscoelastic behavior, this technique would be appropriate for the stress-relaxation calculation; however, for the viscoelastic-creep and mechanical-hysteresis studies, where the system is subjected to specified external stress, an additional dynamic equation is needed.

Pressure-Controlled Dynamic Equation.
The study of viscoelastic creep and mechanical hysteresis involves specifying the external stress, or, in this study, external pressure, and investigating how the (average) internal pressure changes temporally; this requires the dynamic equation that contains external pressure as an independent variable. Such equation, referred to as the pressure-3 Advances in Mathematical Physics controlled dynamic (PCD) equation, was first introduced by Kocher and Provatas [38] expressed as where P ext is the externally applied pressure andL def ′ is the mobility coefficient. This equation can be used in conjunction with Equation (7) to simulate the system under specified P ext . If the external pressure is higher (lower) than the internal pressure, or P ext > P int ð P ext < P int Þ, the value of dṼ/dt will be negative (positive) which indicates that the system will undergo compressive (tensile) deformation.
A slightly modified version of the PCD was recently proposed by Em-Udom and Pisutha-Arnond [37]. The equation is expressed as whereL def is the mobility coefficient. This equation was developed using the principle of thermodynamics of a hydrostatically stressed solid [39,46] and classical irreversible thermodynamic framework [40][41][42]. The comparison between results between Equations (9) and (10) shows that the steady-state behavior obtained from both equations is identical; however, the nonequilibrium behaviors are different. This difference is due to the extra volume term (Ṽ) on the right-hand side of Equation (10). This extra volume term results in a higher deformation rate in tension (due to increasingṼ) than the deformation rate in compression (due to decreasingṼ). Nevertheless, for small deformation, the numerical difference from Equations (9) and (10) is unlikely to change the qualitative interpretation of the results. Despite the fact that we employ Equation (9) in this work, the choice of the PCD equation is immaterial for this study. Hereafter, we will refer to the PFC method that incorporates the PCD equation as the PFC-PCD model. Lastly, the tilde notation will also be omitted for simplicity. More detail on the PFC-PCD model is provided in Appendix A.

Standard Linear Solid Model.
We review the standard linear solid (SLS) model which is a method of modeling viscoelasticity using linear combinations of springs and dampers as shown in Figure 1. The strain response from this model will be compared with that from the PFC-PCD simu-lation in order to extract viscoelastic quantities and properties. The governing differential equation for the SLS model is [20,47] where e is the strain response from the applied stress σ. Equation (11) can be written in simple form as follows: where K = E 1 E 2 /ðηðE 1 + E 2 ÞÞ and τ = η/E 1 . The variable σ ′ is the scaled quantity defined as σ ′ = σ/ðτðE 1 + E 2 ÞÞ. For the creep behavior, we consider the case where the system (initially at zero stress and strain) is subject to constant stress: The corresponding strain is then where e ∞ = σ ′ /K is the deformation at t = ∞ and τ c = 1/K is the relaxation time; the subscript "c" indicates the creep phenomena. The quantity e ∞ is the deformation at steady state while τ c relates to the rate at which the system reaches the steady-state deformation. The plot of the strain response is shown as an example in Figure 2. For the hysteresis behavior, we consider the system induced by harmonic excitation given by

Advances in Mathematical Physics
where σ A ′ is the stress amplitude, ω = 2πf is the angular frequency, and f is the frequency. At steady state, the strain response yields the following solution (see Appendix B): where is the strain amplitude. As a result, the amplitude ratio, A r , is Also, the loss tangent factor, defines the phase angle due to time-delayed response. From the solution, the complex moduli (or dynamic modulus) The quantity G complex ′ is a property of viscoelastic materials and can be written further as where G storage ′ is the storage modulus and G loss ′ is the loss modulus. The quantity G storage ′ measures the stored energy and is an elastic response of material while G loss ′ measures energy dissipation and is a viscous response of a material. The expressions for G storage ′ and G loss ′ can be obtained from rearranging the last equality in Equation (20), yielding The loss tangent, tan ðδÞ, is the ratio of storage and loss modulus or tan ðδÞ = G loss ′ /G storage ′ ; this quantity is a measure of damping in the material.
The behavior of the SLS model, a function of excitation frequency, is shown in Figure 3(a). At low frequencies, tan ðδÞ or δ is low, which indicates that the system behaves in an elastic manner. This is due to fact that the dashpot in Figure 1 has sufficient time to displace, resulting in the stress and strain profiles that are almost in phase. The value of G storage ′ is relatively small compared with that at higher frequencies due to large strain amplitude ( Figure 3(b)). At high frequencies, the system also behaves elastically as seen from small tan ðδÞ; at this condition, the change in stress is too rapid for the dashpot to operate, and thus, the strain response is in phase with the stress profile. Due to small strain amplitude, the value of G storage ′ is relatively large. At intermediate frequencies, a considerable amount of phase lag occurs, as shown schematically in Figure 3(c), showing a more viscous behavior, and the stress-strain profiles can be plotted to exhibit a hysteresis loop as shown in Figure 3(d) where the area inside the loop represents the energy loss. At the frequency where tan ðδÞ is maximum, the system exhibits the highest damping capacity, and in real materials, this condition is important for damping applications.

Method of Simulation
In this section, we describe the setup and method of simulations. In all simulations, we construct a BCC unit cell under hydrostatic pressure ( Figure 4) with a periodic N × N × N grid where N = 16 and the lattice parameter L a = 2 π ffiffi ffi 2 p ; this results in the reference grid spacing, Δr 0 = ffiffi ffi 2 p π/8 in all dimensions, and the reference volume The initial density profile is constructed with the one-mode approximation in Equation (6), where the reference density, ρ 0 , is specified. The profile is then relaxed using the evolution equation (Equation (7)), while fixing Δr = Δr 0 and ρ = ρ 0 , until the equilibrium density profile is reached. At this state, the reference internal pressure, P 0 , can be calculated using Equation (8). Also, the numerical method employed to solve the evolution equation, and also the PCD equation (Equation (10)), is the Fourier pseudospectral method where the Fourier transform is used to calculate spatial derivatives and Figure 4: BCC unit cell under hydrostatic pressure schematic using isosurface construction; σ ii = P ext depicts stress in the fiig direction where all σ ij components, shear stress components fi ≠ jg, are equal to zero.
We note that the grid spacing ΔrðtÞ is identical in all directions due to the assumption of a crystal under hydrostatic stress. During the simulation, the mass conservation is imposed by changing the average density through Also, the mechanical equilibrium is maintained throughout the deformation process; this is accomplished by relaxing the system using the evolution equation at every time step of the time evolution from the PCD equation. The system response is characterized by the strain defined by which is identical in all directions. The profiles of eðtÞ and P ext ðtÞ are then used to quantify the viscoelastic behavior predicted form the PFC-PCD method. Two types of P ext ðtÞ are used in this study. For the viscoelastic creep simulation, the applied pressure is set to where P ext is a constant. The values of P ext are chosen such that the pressure deviation from the reference pressure is positive or negative; the former leads to tensile deformation while the latter results in compressive deformation.
The quantities e ∞ and τ c are numerically extracted from the strain response, eðtÞ, where e ∞ is the values of eðtÞ at large t and τ c is the value of t at eðτ c Þ = e ∞ ð1 − exp ð−1ÞÞ ≈ 0:6321e ∞ (see Equation (14) and Figure 2). For the hysteresis simulation, the applied pressure is set to where P A is the pressure amplitude. The pressure deviation from the reference pressure will then be With this type of pressure function, the system is then allowed to evolve temporally until the steady-state strain response is reached. The strain response profile at steady state will exhibit the sinusoidal behavior similar to Figure 3(a), and we can numerically measure the strain amplitude e A in order to calculate the amplitude ratio

Results and Discussions
4.1. Viscoelastic-Creep Behavior. In this part, we show the viscoelastic-creep behavior exhibited by the PFC-PCD model and demonstrate the dependence of this behavior on Δ P 0 and L def . In Figure 5, we aim to show the strain response from different values of Δ P 0 . The graphs show the strain, initially at zero, increase or decrease gradually until the steady-state values (denoted as e ∞ ) are reached. These delayed responses obtained from the simulations are indicative of the viscoelastic-creep behavior; this means that the PFC model combined with the PCD equation is capable of modeling viscoelastic behavior. Also, depending on the values of Δ P 0 , either tensile or compressive deformations are realized. When Δ P 0 > 0, or P ext < P 0 , tensile deformation is obtained while Δ P 0 < 0, or P ext > P 0 , results in compressive deforma-tion. The magnitudes of e ∞ (steady-state deformation) also scale with the magnitudes of Δ P 0 , as expected from the elastic effect in the PFC model (Equation (1)).
Next, Figure 6 shows how viscoelastic behavior is affected by L def . From the figure, the increase in L def results in eðtÞ approaching e ∞ at a faster rate, or smaller τ c ; this means that the increase in L def leads to the faster strain response or higher degree of elasticity. This result demonstrates that L def can be used to adjust the degree of viscoelasticity ranging from highly viscous to highly elastic. The capability to control the degree of viscoelastic creep can be used to tune the PFC-PCD model to different types of materials.
The degree of viscoelasticity is also affected by ϵ (inverse of temperature), as shown in Figure 7. From the figure, the values of τ c decreases as the values of ϵ increases for the range of P 0 where the solid is stable. Since ϵ is inversely proportional to temperature, this indicates that the system responds faster, or becomes more elastic, as temperature decreases; this prediction agrees with the experimental observation where relaxation time in creep phenomena decreases with temperature.
We note that the results from Figures 6 and 7 assume that L def and ϵ are independent from one another. Nevertheless, since both L def and ϵ affect the viscoelasticity behavior, it is likely that L def and ϵ (or temperature) are related. From the derivation in [37], L def originates from linear phenomenological law in the classical irreversible thermodynamic framework. Therefore, the dependence of L def on temperature could be introduced by considering a more rigorous treatment of this phenomenological law.

Advances in Mathematical Physics
Lastly, we report how the stiffness changes with ϵ and ρ 0 . This is shown in Figures 8(a) and 8(b), where the values of e ∞ are plot as functions of ϵ and ρ 0 , respectively. Since Δ P 0 is the same in all simulations, lower e ∞ values indicate higher stiffness and vice versa. From Figure 8(a), e ∞ decreases as ϵ increases while, from Figure 8(b), e ∞ increases as j ρ 0 j increases. This indicates that the stiffness is higher at lower temperature and lower atomic densities. These results are qualitatively consistent with the theoretical framework [13,48] and experimental observation [49][50][51]. We note that stiffness is not a viscoelastic property and the consistency of the stiffness variation with ϵ and ρ 0 originates from the PFC free energy (Equation (1)), not the PCD equation. Nevertheless, these stiffness results demonstrate the application of the PCD equation since it allows convenient calculation of e ∞ from specified values of Δ P 0 . The results indicate that the PFC-PCD model does not exhibit elastic behavior at very frequency. This suggests that even though the PFC-PCD model is capable of exhibiting viscoelastic response (or phase lag) from the oscillating pressure, the applicability of these phenomena is still limited to the low-frequency excitation.
Another comparison between the PFC-PCD and the SLS model can be shown in Figure 10 where the amplitude ratios are shown. Both the PFC-PCD and the SLS model show reducing amplitudes with increasing frequencies and both amplitude ratios exhibit qualitatively similar functional forms. However, the amplitude ratio from the PFC-PCD model is vanishingly small at high frequencies while the amplitude ratio from the SLS model remains at finite values. The values of the amplitude ratio can also be seen from the hysteresis loop in Figure 11 where the results in Figure 9 are plotted on stress/pressure-strain axes. As the frequency increases, the hysteresis loop rotates counter-clockwise, indicating smaller strain amplitude (or amplitude ratio). At f = 100, the hysteresis loop from the PFC-PCD model almost forms a vertical line due to an essentially zero amplitude ratio while the hysteresis loop from the SLS model remains tilted away from the y-axis. This extremely small amplitude ratio indicates that the system produces almost zero strain from the input pressure, which leads to unnaturally large values of the moduli (G′ storage and G′ loss ). Therefore, the application of the PFC-PCD method to model dynamic viscoelastic behavior should be limited to low-frequency excitation. Nevertheless, at the low frequency, the prediction of G ′ storage , G′ loss , and tan ðδÞ is in good agreement with the results from the SLS model as shown in Figure 12.
To further investigate into the difference between the results from the PFC-PCD and those from the SLS model, we expand the PCD equation in terms of strain and perform Taylor's expansion to the first order with respect to the strain; the details of the derivation are shown in Appendix C. Compared with the SLS equation, the transformed PCD equation lacks the term with the time derivative of stress. This analysis indicates that the PCD equation is more similar to the Kelvin-Voigt model than the SLS model. The Kelvin-Voigt model consists of a spring and a dashpot in parallel and is capable of reproducing the creep phenomena, but not the stress relaxation. This similarity between the PCD and the Kelvin-Voigt models explains why the PFC-PCD model can reproduce creep phenomena, but not the dynamic behavior similar to the SLS model. We note that the lack of similarity to the SLS model not only applies to our recently proposed PCD equation [36] but also applies to the PCD equation proposed by Kocher and Provatas [38] since the functional Similar to the result from the viscoelastic-creep study, the PCD parameter L def can be used to control the degree of viscoelasticity. Figures 13(a) and 13(b) show the strain response with different values of L def at a fixed pressureoscillation frequency. The hysteresis loops from these results are also shown in Figure 14. The results show that decreasing the value of L def leads to higher phase lag ( Figure 13) and higher damping capacity ( Figure 14). This indicates that lowering the L def value leads to the system having a more viscous response. The opposite is also true where increasing L def leads to a more elastic response from the system.
The influence of ϵ on G′ storage , G′ loss , and tan ðδÞ is shown in Figure 15. Figure 15(a) shows that an increase in ϵ (decrease in temperature) leads to an increase in G ′ storage and a decrease in G ′ loss . As tan ðδÞ = G ′ loss /G ′ storage , the values of tan ðδÞ decrease with increasing ϵ, as shown in Figure 15(b). The results show that as temperature decreases, the system behaves in a less viscous manner and with reducing damping capacity. These trends qualitatively agree with polymeric materials below glass transition [52] and metals such as aluminum alloy [53].
It should be noted that the particular form of the PFC free energy (Equation (1)) is applicable to crystalline materials such as metals [45]; therefore, it might seem that the results from this work only pertain to these types of materials. Nevertheless, since the viscoelastic behavior originates from the PCD equation, not the PFC free energy, an alternative form of PFC free energy (Equation (1)) can be used and the modified PFC-PCD model should still exhibit similar viscoelastic behavior as well as similar dependence of viscoelastic properties on the PFC-PCD model parameters. The capability to generalize the PFC free energy allows modeling of viscoelastic behavior in different types of materials such as polymeric materials where viscoelasticity is much more pronounced than that in metals.

Conclusion
In this study, we aim to show the capability and the limitation of the PFC-PCD model which was first employed to

12
Advances in Mathematical Physics investigate the viscoelastic-creep and mechanical-hysteresis behaviors in a BCC (3D case) unit cell under hydrostatic pressure. To achieve our goal, we implement two types of pressure profiles in our analysis: constant pressure input for viscoelastic-creep study and sinusoidal pressure oscillation for mechanical-hysteresis study where all PFC parameters, temperature and atomic density parameters, are established in a solid stable region. The conclusions for each analysis are summarized as below: 1. In the viscoelastic-creep study, we first investigate the influence of two PCD parameters: magnitude of pressure input and mobility coefficient on viscoelastic-creep behavior. Regarding the results, we find the magnitude of pressure input scales with that of deformation at steady-state time both tensile and compressive loads. The mobility coefficient can be used to control the degree of viscosity of deformation in that higher mobility results in higher elasticity as well as lower mobility leads to higher viscosity. This evidence suggests that we can control both the magnitude of deformation and the degree of viscoelasticity through PCD parameters. Second, we study the impact of PFC parameters: temperature on relaxation time and the magnitude of deformation at steady-state condition as well as atomic density parameter on the magnitude of deformation at steady-state condition.
Regarding the results, we find that the system response is faster and has less relaxation time, less deformation, and more stiffness, at lower temperature which is generally consistent with theoretical framework and experimental observation. Also, the system exhibits larger deformation and less stiffness, at higher atomic density that shows an agreement with previous work 2. In dynamic behavior study, we first focus on studying the impact of excitation frequency on strain response, amplitude ratio, mechanical hysteresis, and complex moduli. The results predicted by the PFC-PCD model indicate that the system displays more phase lag between sinusoidal pressure oscillation and strain response with reducing amplitude ratio at increased frequency which qualitatively agrees with that of the SLS model. However, at high frequency range, the results exhibited by both models are inconsistent. This finding can also be observed in the mechanicalhysteresis loop. The hysteresis loop exhibited by PFC-PCD becomes a vertical line which indicates the system produce almost zero strain. This finding produces unusually large values of complex moduli unlike the SLS model prediction. Furthermore, we show that the mobility coefficient can also be used to control the degree of viscoelasticity in dynamic behavior under sinusoidal pressure oscillation. The impact of temperature on complex moduli predicted by the PFC-PCD model is also in good agreement with experimental observation Finally, we show that the PCD model is more similar to the Kelvin-Voigt model than the SLS model regarding Taylor's expansion in Appendix C. This deviation does not seem to cause any difference in viscoelastic-creep behavior, but the mechanical-hysteresis behavior tends to differ from the SLS model at high frequency range. However, the prediction of the model agrees well with the SLS model at the lowfrequency range. Therefore, the application of the PFC-PCD method to predict dynamic viscoelastic behavior should be limited in some certain frequency region. This analysis provides the possibility of enhancing the PFC-PCD model by modifying the PCD equation in a more proper way which might be addressed in future improvement.

Advances in Mathematical Physics
The PFC-PCD equation consists of two governing equations: Cahn-Hilliard-(CH-) type equation and the pressurecontrolled dynamic (PCD) equation. The first equation is involved in dissipative dynamics and a mass conservation equation [2] formulated in partial differential equation (PDE) with four independent variables ðt,r 1 ,r 2  In order to solve Equation (A.10) numerically, there are two options which can be employed: finite difference method (FD) and Fourier spectral method (FSM). In this study, not only is the periodic boundary condition governed but also the equation itself has a very high derivative order in spatial coordinates which demands a lot of numerical expense if the FD scheme is used. Therefore, the FSM method is more appropriate than FD in terms of numerical accuracy and numerical expense to solve Equation (A.10). To employ the FSM scheme, we begin with Equation (A.9). If we apply the Fourier transform both sides, we obtain ðA:11Þ Then, the 6OD-NPDE is literally reduced to an empirical ordinary differential equation which depends on time constraint only, where b ρ ≡ b ρðK,tÞ is the Fourier transform ofρ.ρ∧ 3 ≡ρ∧ 3 ðK,tÞ is the Fourier transform ofp 3 which can be described as follows: where d is dimensionality of the problem. This variable will be used to calculate P int in the pressure-controlled dynamic equation further. Next, the second important equation is the PCD equation which was developed regarding the thermodynamics of the hydrostatically stressed crystal solid [39] and classical irreversible thermodynamic frameworks [40][41][42]. The advantage of this equation is that it allows creating the deformation in the PFC model by specifying the external pressure as an input variable unlike the conventional PFC simulations where the volume, or grid spacing, is an input variable [43]. This equation is given as follows:

B. Standard Linear Solid Model
In this appendix, we aim to derive the expressions for the complex modulus and strain response from the SLS model under harmonic excitation. First, we derive complex modulus under dynamic excitation. Consider a more general form of the stress (scaled quantity): where σ A ′ is the amplitude, ω is the angular frequency and σ ′ * A is the complex amplitude defined by σ ′ * A ≡ σ A ′ exp ðjδÞ = σ′ * A cos ðδÞ + jσ′ * A sin ðδÞ. At steady state, the strain response will also exhibit a harmonic function with the same frequency; therefore, one can postulate the functional form of the strain as an input function [20]:  It should be noted that e A = je * A j.

C. Dynamic Similarity between PCD Equation and SLS Model
In this appendix, we demonstrate how dynamic similarity between the PCD Equation (10) and the SLS model Equation (12) can be established using Taylor's expansion [55] at infinitesimal deformation around the reference state. First, we begin with PCD equation where P int = Ð V P int dV/V and P ext represents the average internal pressure and external applied pressure with VðtÞ = LðtÞ 3 where L is the crystal length parameter. If we denote If we assume infinitesimal deformation, small strain can be written as eðtÞ = e 0 + Δ e where e 0 = 0 and Δe = eðtÞ − e 0 . Therefore, Equation (C.3) can be written in expansion form [55]: where φðe 0 Þ represents Ð V PdV at e = e 0 and represents a higher order derivative term as a function of Δe 3 . In case of infinitesimal deformation, the higher order term, Δe n where n ≥ 2, can be truncated, then Equation (C.4) is reduced to φ e ð Þ ≈ φ e0 ð Þ + ψ e=e 0 Δe, ðC:5Þ where ψ e=e 0 = dφ/deje = e 0 . Also, P ext can be written as P ext = P 0 − Δ P ext , ðC:6Þ where P 0 represents average reference pressure calculated under an undeformed state P 0 = φ e 0 ð Þ L 3 0 , V 0 = L 3 0 ,

ðC:7Þ
and Δ P ext represents the deviation part of P ext from P 0 . Since V ≡ VðLðeÞÞ is a composite function, therefore, the time derivative regarding Equation (C.1) can be expressed by the chain rule  Regarding binomial expansion, we know that ð1 + eÞ −2 ≈ 1 − 2e for the first term approximation if 0 ≤ e<<1. Equation (C.11) becomes which obviously has structure like a Kelvin-Voigt model: where K ≡ L P ð3 P 0 − ψ e=e 0 /L 3 0 Þ/3 and σ ′ ≡ L P Δ P ext /3. Equation (C.15) is slightly different from the SLS model Equation (12) in that the τdσ′/dt term is introduced on the right-hand side in the SLS model. This deviation does not cause the difference in viscoelastic-creep behavior produced by both the PCD model and the SLS model solution under constant pressure. However, in case of the mechanical hysteresis, although this behavior produced by the PCD equation is deviated with that of the SLS model and behaves like the Kelvin-Voigt model at an elevated frequency range, it consistently predicts similar results with the SLS model at low frequency. This can be described by the solution obtained from both models under the sinusoidal load, σ ′ ðtÞ = σ ′ A sin ðωtÞ. According to Equation (B.9), the strain solution produced by the SLS model is eðtÞ = σ A ′ sin ðωt − ϕÞ/ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K 2 + ω 2 p + τσ ′ A ω cos ðωt − ϕÞ/ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K 2 + ω 2 p . However, if 0 ≤ ω<<1, the solution will converge to eðtÞ = σ A ′ sin ðωt − ϕÞ/ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi K 2 + ω 2 p which is almost identical to that of a simplified model of the PCD equation.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.