Kinetics of orientational phase ordering near line defects in crystals

KINETICS OF ORIENTATIONAL PHASE ORDERING

NEAR LINE DEFECTS IN CRYSTALS

C. Bjerk´ena_{, A. R. Massih}a,b

a_{Fakulteten f¨or teknik och samh¨alle, Malm¨o h¨ogskola, SE-205 06 Malm¨o, Sweden} b_{Quantum Technologies AB, Uppsala Science Park, SE-751 83 Uppsala, Sweden}

Keywords: Phase transitions, Kinetics, Defects, Ginzburg-Landau theory Abstract

General properties of directed ordering near line defects in elastic crystals undergo-ing phase transition are studied usundergo-ing the two-component time-dependent Ginzburg-Landau equation. Upon quenching the system below its transition point, the temporal evolution of the order parameter in the vicinity of the defect is evaluated. The de-velopment of vortices is explored and their interaction with the structural defect is examined. Finally, phase transitions in improper ferroelectrics in the context of the model are discussed.

Introduction

Oriented ordering of second-phase occurs in many metallurgical systems. For exam-ple, in transition metals, such as titanium and zirconium, hydride precipitates assume different orientations depending on the material texture and the direction of the ap-plied stress field [1]. In phase-field modeling, for such a system, one assumes that the order parameter is a vector field whose components describe differently oriented variants of the precipitate phase. Besides phase precipitation, vector order param-eter field has been used to recount a variety of phenomena in nature [2]. Here, we study some basic properties of phase ordering near line defects in elastic crystals un-dergoing phase transition using the two-component time-dependent Ginzburg-Landau equation, which describes two different low symmetry phases. The Landau potential comprises the elastic properties, the singularity that characterizes the defect and a bi-quadratic anisotropy term producing different symmetries. The study is an extension of [3] which considered the order parameter as a scalar field variable.

Theoretical Setup

We set up the Ginzburg-Landau equation for modeling situations, where in a crys-talline solid embracing a line defect, different low-symmetry phases would emerge upon quenching the system from its high-symmetry phase below its transition tem-perature. This system is described by a two-component order parameter η = (η1, η2) obeying the time-dependent Ginzburg-Landau (TDGL) equation:

∂ηi

∂t = −Γ δF δηi

, i = 1, 2, (1)

1_{Presented at PTM 2015, June 28-July 3, 2015, Whistler, BC, Canada.}

where Γ is a kinetic coefficient and F is the total energy of the system F [η] = Z hg 2 2 X i=1 (∇ηi)2+ V(η) i ddx. (2)

Here, g is taken to be a positive constant, d is the spatial dimensionality of the system and V(η) is the two-component Landau potential expressed as

V(η) = α 2|η| 2₊ β 4|η| 4₊λ 6|η| 6₊ 1 2 γ1+ γ2|η| 2
_η2 1η22, (3) where |η| =pη2

1 + η22 and α, β, λ, γ1, γ2 are phenomenological parameters specified later. In particular, for the problem under consideration, α is a space-dependent quantity including the strain field produced by a defect. Note that the cross term η2

1η22 in (3) breaks the SO(2) symmetry of the potential. The equations for the space-time variation of the order parameter components are now obtained by inserting (2)-(3) into (1) 1 Γ ∂ηi ∂t = g∇ 2_η i− ∂V(η) ∂ηi , i = 1, 2 (4) where ∂V(η) ∂ηi =α + β|η|2+ λ|η|4+ γ1ηj2+ γ2(η2iη 2 j+ η 2 j|η| 2_)η i, i 6= j(= 1, 2). (5) As alluded above, in the presence of a defect producing a field in an elastic lattice, e.g. an edge dislocation or a domain wall, the coefficient α becomes space-dependent and the coefficient β would depend on the elasticity parameters [3, 4, 5]. More precisely, we may write: α(x) = a τ + u(x)
, β = β0− 2ς2/Λ0, and τ = (T − Tc)/Tc, where Tc is the global Curie temperature, a taken as a positive constant, ς the striation coefficient, which accounts for the strength of the interaction between the elastic field and the order parameter, Λ0 = K + (2 − 2/d)M , K the bulk modulus, M the shear modulus, and u(x) denotes the strain field of the defect. For a defect free crystal α = aτ , while for a rigid crystal β = β0. For an edge dislocation embedded in a crystal the strain field, in polar coordinates x = r(cos θ, sin θ), is u(r, θ) = B sin θ/r, where B = _2π(1−ν)b(1−2ν), b is the magnitude of the Burgers vector and ν is the Poisson ratio [4]. In case of a crystal embracing a coherent twin boundary, a strain field in the form u(x) = A cosh−2(x/x0) contributes to the energy of the system, where A is a material constant equivalent to B and x0( 1) is the thickness of a thin twin boundary; thus making u(x) short-range [6].

Steady State, Equilibrium and Phase Diagram

In a steady state the right-hand side of (1) is zero, hence the order parameter becomes time-independent. Thus (4) and (5) give

g∇2ηi =E + au(x) + β|η|2+ λ|η|4+ γ(ηi, ηj)ηi, i 6= j(= 1, 2), (6) where the boundary conditions are the absence of singularity at |x| = 0 and a finite limit as |x| → ∞. Here, we put: γ(ηi, ηj) ≡ γ1η2j + γ2(η2iηj2+ η2j|η|2) and E ≡ aτ . Now in line with arguments in [4] there exists a certain E = E0 < 0, corresponding to T0 > Tc, such that (6) has a nontrivial solution for E > E0(⇒ T < T0). The

value of E0(≡ a(T0− Tc)/Tc) may be determined from the solution of the linearized equations, which provide the state eigenfunctions of the problem. A bound-state solution η(x) exists if F [η] < 0. For E > E0, if we substitute the eigenfunction η = qψE0(x) of the linearized equation (6) in (2), corresponding to its eigenvalue

E = E0 where q ≡ qE0 < 1, we find F [qψE0] < 0. For a two-component order

parameter, we simply write ηi(x) ≈ qiψE0(x), i = 1, 2. Inserting this ηi(x) in (2),

using g∇2_ψ

E ≈ [E + au(x)]ψE, and integrating over the volume, the energy of the system can be written as

F (q1, q2) = F0+ α∗ 2 |q| 2₊ β 4|q| 4₊λ 6|q| 6₊1 2γ(q)q 2 1q22, (7)

where F0 is the energy of the high-temperature phase, α∗ = E − E0, |q| =pq12+ q22, γ(q) = γ1 + γ2|q|2, and we have absorbed the constants resulting from the volume integrations in the coefficients of F .

The possible stable phases associated with (7) or (3), α ≡ α∗, are as follows: The high symmetry phase I (η1 = η2 = 0), α > 0; the low symmetry phases II (η1 6= 0, η2 = 0 or v.v.), α < 0, γ1 > 0 and III (η1 = η2 6= 0), α < 0, γ1 < 0; and possibly phase IV (η1 6= η2 6= 0), α < 0, if γ1 = −γ2. Figure 1 illustrates the three stable phases of the considered system. Minimization of F in (7) or V in (3) gives the phase boundaries (stability lines), namely, α = 0, α = B2_{/4Λ, and} α = B(γ1/γ2) − Λ(γ1/γ2)2 if γ1 < 0 ∧ γ2 > 0 ∨ v.v., where B = β, Λ = λ for phase II or B = β + γ1/2, Λ = λ + 3γ2/4 for phase III.

In the phase diagrams (Fig. 1) the (α∗, γ1)-plane is constructed for β > 0 and β < 0, separately, in which the stability limits for phase I (dotted line), phase II (dash-dot line) and phase III (dashed) are shown. The diagrams show the stability limits and the transition borders. They depict the regions of phase stability overlap, where F has both global and local minima; thus stable and metastable phases may coexist. The second-order phase transition only takes place from I→II or III for β > 0. When there is a range of phase coexistence, a first-order phase transition would occur, as indicated by solid lines in Fig. 1. The upper border (solid) lines are found by the condition F = F0, which gives α∗ = 3B2/16Λ. The II ↔ III transition line is obtained by setting FII= FIII, where the resulting equation is solved numerically. Recall that for β < 0 (λ > 0, γ2 > 0) phase transitions are first-order, the origin α∗ = β = γ1 = γ2 = 0 corresponds to the tricritical point, and note that the transition borders depicted in Fig. 1 meet at a single point in the (α∗, γ1)-plane, viz. α∗ = 3β2_{/16λ, γ}

1 = −2β 1 −pΛ/λ
.

Numerical Computations Method

We have used the open source partial differential equation solver package FiPy [7], based on a finite volume method, for our numerical computations. The spatiotemporal evolvement of the two-component order parameter η = (η1, η2) in 2D is evaluated by solving numerically equations (4)-(5). An edge dislocation is placed at the center of a square cell, (x, y) = (0, 0), which has its slip direction along the the positive x-direction. The side length of the cell L = 4r0, where r0 ≡ B. The cell comprises at least 200 × 200 equally sized elements. At the outer boundary of the cell we impose the Neumann condition n · ∇η = 0, where n is a unit vector perpendicular to the

−3 −2 −1 0 −2 −1 0 γ₁ / β α *

I

II

III

I

II

III

Figure 1: Phase diagrams relating to the thermodynamic potential (7) for (a) β > 0, and (b) β < 0, showing stability limits and transition borders; λ > 0 ∧ γ2 > 0.

boundary. We take an initial value of η to be a tiny random number within ±5 × 10−10 and set our reference time step to ∆tref = 0.01, which is also our actual computation time increment. A numerical simulation comprises 104 _{time steps.}

Results of our computations for the evolution of order parameter can be presented either in a Cartesian form η = (η1, η2) or in a polar form η = ρ(cos ϕ, sin ϕ). The evaluate the evolution of the phases, we set ϕ = nπ/2 for phase II and ϕ = (2n+1)π/4 for phase III, where n is an integer. Phase I, which is assumed to predominate the entire system prior to the insertion of the dislocation, is defined by setting ρ = 0. Results

Evolution of order parameter. We investigate the evolution of η by choosing different combinations of the coefficients in (3). We present the results for β = γ1 = −1, aτ = λ = γ2 = 1 and Γ = g = 1. We scale down η by η0 = (|aτ |/|β|)1/2, i.e. η/η0, to get a dimensionless quantity, where η0 has the same order of magnitude as the maximum of η obtained by our computations. Figure 2 shows our results of the spatial distribution of the order parameter components η1 (upper row) and η2 (lower row) at two different times t = 103_∆t

ref and 104∆tref, respectively, in the form of contour plots. It is seen that close to the dislocation, η1 reaches its maximum value and the area where η1 deviates from zero markedly spreads out until it ceases at ≈ 9000∆tref. A different behavior is seen for η2 where kidney-shaped contour levels evolve and their peak keeps rising, but never reaches the level of η1, besides it is further located from the position of the dislocation. This means that close to the dislocation η1 6= 0, while η2 ≈ 0. Further away, both η components are non-zero, although not equal. This may indicate that two different low-symmetry phases are forming for this particular parameter set (β = γ1 = −1).

The associating plots of (η1, η2) versus r/r0 in various times are depicted in Fig. 3 which are for times from t = 2 × 102∆tref to 104∆tref. As can be seen, first both η1 and η2 rise in the vicinity of the dislocation, but η1 grows faster then swiftly forms a sharp top with a peak value at a position that remains fairly constant thereafter, while its width broadens. The evolution of η2 is different, i.e. the initial top will grow

somewhat until a peak value of η2/η0 ≈ 0.5 is reached. Subsequently, the top reshapes to form a kidney shape figure with a maximum in the (x, y)-plane at y/r0 ≈ −1, cf. Fig. 2(d). −1 0 1 −2 −1 0 x / r₀ y / r 0 (a) η1/η0 at t = 1000∆tref −1 0 1 −2 −1 0 x / r₀ y / r 0 0.5 1 1.5 2 (b) η1/η0at t = 10000∆tref −1 0 1 −2 −1 0 x / r₀ y / r 0 (c) η2/η0at t = 1000∆tref −1 0 1 −2 −1 0 x / r₀ y / r 0 0.1 0.2 0.3 0.4 0.5 (d) η2/η0at t = 10000∆tref

Figure 2: Spatial distribution of the components of η at two different times for aτ = λ = γ2 = 1 and β = γ1 = −1; n.b. the color-bar scales are different for the two components.

Vortex patterns. A basic attribute of the 2D two-component Ginzburg-Landau the-ory, or the 2D XY model, is the presence of vortices (antivortices) below a transition temperature. To realize or visualize the evolution of the vortex patterns in the system under consideration, we compute the distribution of the order parameter phase angle ϕ = arctan(η2/η1). In these computations, we use a denser mesh, 400 × 400, and a smaller time step, 0.1∆tref, than previously to clearly resolve the details in variations of ϕ(x, y) at different times. The results of some of our computations are displayed in Fig. 4 as contours of the phase angle ϕ, or rather sin2(2ϕ), to conveniently produce them as greyscale images. Figs. 4(a)-(c) show the vortex/antivortex patterns at three different stages of the evolution for a cell with no dislocation, whereas Figs. 4(d)-(f) display the corresponding results with the presence of a dislocation. In these compu-tations, we have used β = 1, γ1 = −1, γ2 = 1 and assumed aτ = λ = Γ = g = 1. A pattern of vortices emerges from the initial perturbation of the order parameter and then coarsens with time due to vortex-antivortex coalescence or annihilation. By in-troducing a dislocation in the system, the vortex pattern is affected by the dislocation field such that an area develops and expands around the dislocation where the vor-tices are annihilated efficiently. A more detailed analysis of vortex pattern evolution

0 1 2 0 0.5 1 1.5 2 r / r0 η 1 / η 0 0 1 2 0 0.5 1 1.5 2 r / r₀ η 2 / η 0 t / ∆ t_ref= . . . 200 400 600

Figure 3: Evolution and spatial distribution of the order parameter components (η1, η2) around a dislocation (r = 0) in various times in sequence {n × 102∆tref}; n = 2, 4, . . . , 100. Here θ = −π/2. −2 −1 0 1 2 −2 −1 0 1 2 (a) t = 50∆tref −2 −1 0 1 2 −2 −1 0 1 2 (b) t = 100∆tref −2 −1 0 1 2 −2 −1 0 1 2 (c) t = 500∆tref −2 −1 0 1 2 −2 −1 0 1 2 (d) t = 50∆tref −2 −1 0 1 2 −2 −1 0 1 2 (e) t = 100∆tref −2 −1 0 1 2 −2 −1 0 1 2 (f) t = 500∆tref

Figure 4: Snapshots of vortex patterns in the (x, y)-plane, where Schlieren patterns with intensity sin2(2ϕ), ϕ = arctan(η2/η1), develop in various times for γ1 = −1 and γ2 = 1: (a)-(c) No dislocation, (d)-(f) with dislocation. Greyscale: Black 0 − 0.25; dark-grey 0.25 − 0.5; light-grey 0.5 − 0.75; white 0.75 − 1.

comprising the effect of noncubic versus the SO(2) symmetric potential (γ1 = γ2 = 0) is discussed in [8].

Discussion

In Fig. 1, we have depicted diagrams separating the aforementioned phases for the system. It may be worth mentioning that according to [9], phase I has the point group symmetry C4v (tetragonal), phases II and III have point group symmetry C2 (monoclinic) and any possible domains of phases II and III correspond to a domain wall whose symmetry can be either C2 or C1 (triclinic). These symmetry classes, a manifest of the thermodynamic potential (7), may be referred to as “noncubic” symmetry [10]. Physical realizations of the aforementioned phases may be found, e.g. in improper ferroelectrics [11], or in the second-phase precipitation phenomenon. Phase I corresponds to solid solution and II and III may represent different orientations or structures of the precipitates.

The spatial variation of the order parameter components at several times shown in Fig. 3 indicates a kind of self-similarity. Is there a scaling law that would collapse these curves into a single curve for each component? More precisely, could one find a time-dependent characteristic length ξ(t) which scales the distance such that it would lead to a dynamic scaling law [12]? This possibility has not been explored here.

It is known that with the considered model but with SO(2) symmetric potential and no dislocation (2D XY model) if the system is quenched instantaneously below the Berezinskii-Kosterlitz-Thouless temperature [13], vortex-antivortex pairs nucleate, grow and coalesce, thereby reducing the vortex density with time [14]. Indeed, our computations indicate (cf. Fig. 4) that for the case of noncubic Landau potential (γ1 = −1, γ2 = 1), with or without dislocations, the vortex number density evolves with time roughly as % ∼ 1/t [8]. Hence the aforementioned characteristic length evolves as ξ(t) ∼ t1/2. Furthermore, we have found that the diameter of the vortex free area around the dislocation, Figs. 4(d)-(f), increases linearly with time in early times [8].

A direct application of the present model is to the case of the domain structure of improper ferroelectrics such as Gd2(MoO4)3 or GMO [11]. In GMO a phase transition of displacing type occurs at about 430 K upon which the paraelectric tetragonal phase transforms to the ferroelectric orthorhombic (ordered) phase, manifesting the change in the space group symmetry, D3

2d → C2v3 . The phase transition is accompanied by doubling of the unit cell plus an appearance of a two-component order parameter (η1, η2) describing the amplitudes of the lattice distortions. The homogeneous phase of GMO is quadruply degenerate and possess four ground states corresponding to four points in the order parameter plane (η1, η2). The Landau potential used for GMO is isomorphous to (3). Another example is the structural transition in lead zirconate PbZrO3 [15] with a Landau potential isomorphous to (3). In [8], we have considered the case of phase transition in GMO in the presence of a thin twin-boundary contributing a δ-function potential (≈ cosh−2x) to the free energy. An increase in the transition temperature in the vicinity of a twin-boundary is ∆Tc ∝ A2/ag, which results in up to about 5 K for a range of selected model parameters. A cluster of twin-boundaries is expected to yield a higher temperature increase.

References

[1] A. A. Motta and L-Q. Chen. Hydride formation in zirconium alloys. JOM, 64:1403–1408, 2012.

[2] P. C. Hohenberg and A. P. Krekhov. An introduction to the Ginzburg-Landau theory of phase transitions and nonequilibrium patterns. arXiv:1410.7285v2, 2015.

[3] C. Bjerk´en and A. R. Massih. Phase ordering kinetics of second-phase formation near an edge dislocation. Philosophical Magazine, 94:569–593, 2014.

[4] V. M. Nabutovskii and V. Ya Shapiro. Superconducting filament near a disloca-tion. Soviet Physics JETP, 48:480–485, 1978.

[5] A. A. Boulbitch and P. Tol´edano. Phase nucleation of elastic defects in crystals undergoing phase transition. Physical Review Letters, 81:838–841, 1998.

[6] A. A. Bulbich et al. Nucleation at twin boundaries in crystals undergoing phase transitions near singularities in phase diagrams. Soviet Physics JETP, 71:619– 625, 1990.

[7] J. E. Guyer, D. Wheeler, and J. A. Warren. FiPy: Partial differential equations with Python. Computer Science & Engineering, 11:6–15, 2009.

[8] C. Bjerk´en and A. R. Massih. Oriented ordering near line defects. In preparation. [9] A. A. Bulbich and Y. M. Gufan. Inevitable symmetry lowering in a domain wall

near a reordering phase transition. Soviet Physics JETP, 67:1153–1157, 1988. [10] N. W. Ashcroft and N. D. Mermin. Solid State Physics. Harcourt, New York,

1976.

[11] V. Dvorak. Improper ferroelectrics. Ferroelectrics, 7:1–9, 1974.

[12] F. Rojas and A. D. Rutenberg. Dynamical scaling: The two-dimensional XY model following a quench. Physical Review E, 60:212–221, 1999.

[13] J. V. Jos´e. 40 Years of Berezinskii-Kosterlitz-Thouless Theory. World Scientific, Singapore, 2013.

[14] B. Yurke, A. N. Pargellis, T. Kovacs, and D. A. Huse. Coarsening dynamics of the XY model. Physical Review E, 47(1):1525–1530, 1992.

[15] A. K. Tagantsev et al. The origin of antiferroelectricity in PbZrO3. Nature Commun., 4:2229, 2013.