text:
# Self-driven lattice-model Monte Carlo simulations of alloy thermodynamic properties and phase diagrams

## 1 Introduction

Monte Carlo (MC) simulations  of lattice models are a widely used method to compute thermodynamic properties of substitutional alloys. A lattice model MC simulation only requires, as an input, a so-called *cluster expansion* , which defines the energetics of the system by specifying the energy associated with any atomic arrangement on a given lattice. Cluster expansions are typically constructed from a fit to the energies of a set of structures calculated, for instance, from first-principles density-functional-theory calculations , thus enabling the prediction of thermodynamic properties even in the absence of experimental data.

In its simplest form, a cluster expansion specifies the energy difference between nearest-neighbor chemical bonds joining identical and distinct species in a binary alloy ($`E_{AB}\left(E_{AA}+E_{BB}\right)/2`$) and predicts compound formation energies from the number of each type of nearest-neighbor bond. However, this formalism can be extended  to allow for arbitrary interaction ranges and multibody interactions so that the energetics of the alloy can be modeled with an arbitrary accuracy by including a sufficient number of such interactions. It has been demonstrated that the energetic contributions arising from displacements of the atoms away from their ideal lattice sites  can be accounted for within this framework. Moreover, by allowing temperature-dependent interactions, entropic contributions arising from vibrational and electronic excitations can be included as well (see  for a review and  for recent examples). Since long-range interactions are often needed to properly describe the energetics of an alloy system, MC simulations have become a preferred method to compute thermodynamics properties from cluster expansions, because it would be difficult to account for such long-range interactions within accurate mean-field methods, such as the Cluster Variation Method (CVM)  while efficient algorithms to include even infinite-range interactions in MC simulations have been devised . MC simulations are also conceptually simple to understand, straightforward implement, and can be easily adapted to compute a wide variety of properties.

Given these attractive features, it is not surprising that MC simulations of lattice models with cluster expansions determined from first-principles calculations have been employed in a variety of context for metallic, semi-conductor and ceramic systems, including the computation of: composition-temperature phase diagrams, thermodynamic properties of stable and metastable phases, short-range order in solid solutions, thermodynamic properties of planar defects (including surfaces or antiphase and interphase boundaries), and the morphology of precipitate microstructures .

With the steadily increasing availability of computational resources, MC simulations of lattice models are likely to become even more prevalent in computational material science modeling. Despite the technological advances in computational hardware, the human time needed to select the appropriate simulation parameters and to analyze the data has remained essentially constant over the years. We are now in a situation where the main limitation to the inclusion of atomistic MC simulations in standard thermodynamic software toolkits is not the lack of computational resources, but rather the difficulty of controlling a MC simulation in order to obtain the desired quantities.

To address this problem, we have devised a variety of high-level algorithms that serve as an interface between the user and a traditional MC code. The user specifies the goals sought in a high-level form that our algorithm converts into elementary tasks to be performed by a standard MC code. For instance, our algorithms permit the determination of the free energy of a phase over all of its region of stability within a specified accuracy, without requiring any user intervention during the calculations. The proposed algorithms have been implemented in an easy-to-use software package .

In an effort to make this article self-contained, we first recall the basic theoretical results underlying semi-grand-canonical MC simulations, which are especially suited for determination of thermodynamic properties of alloys. We then proceed to describe the algorithms that enable the automated calculation of thermodynamic quantities from MC simulations. These algorithms address the following three questions. (i) How to determine when a MC simulation has run for a sufficiently long of time to provide the desired quantities within a given target precision? (ii) How to detect phase transitions? (iii) How to efficiently determine the composition-temperature phase boundaries of a phase diagram? We then provide a series of examples of applications that illustrate the features of our algorithms.

## 2 Semi-Grand-Canonical Monte Carlo

A very convenient way to obtain thermodynamic information regarding an alloy system is to perform MC simulations which sample a semi-grand-canonical ensemble, also known as the transmutation ensemble. In this ensemble, the energy and concentration of an alloy with a fixed total number of atoms ($`N`$) are allowed to fluctuate while temperature and chemical potentials are externally imposed. Although the results of this paper can be extended to general multicomponent systems, we are considering the case of binary alloys, for the benefit of notational simplicity. In an $`AB`$ alloy, without loss of generality, only the composition $`x`$ of element $`A`$ and the *difference* in chemical potential between the two species $`\mu =\mu _A\mu _B`$ needs to be specified. For conciseness, we simply refer to the quantity $`\mu `$ as the “chemical potential” and $`x`$ as the concentration.

The natural thermodynamic potential $`\varphi `$ (expressed per atom) associated with the semi-grand-canonical ensemble can be defined in terms of the partition function of the system as follows:

$$\varphi (\beta ,\mu )=\frac{1}{\beta N}\mathrm{ln}\left(\underset{i}{}\mathrm{exp}\left(\beta N\left(E_i\mu x_i\right)\right)\right)$$
(1)
where $`E_i`$ and $`x_i`$ are, respectively, the internal energy (per atom) and the concentration of state $`i`$, while $`\beta =1/\left(k_BT\right)`$ is the reciprocal of the temperature $`T`$ and $`k_B`$ is Boltzmann’s constant. The thermodynamic potential $`\varphi `$ is related to the Helmholtz free energy $`F`$ through $`\varphi =F\mu x`$.

The potential $`\varphi `$ can also be defined by the following total differential:

$$d\left(\beta \varphi \right)=\left(E\mu x\right)d\beta \beta xd\mu .$$
(2)
where $`E`$ and $`x`$ are the system’s average internal energy (per atom) and concentration. Equation (2) enables the calculation of the thermodynamic function $`\varphi (\beta ,\mu )`$ through thermodynamic integration:

$$\varphi (\beta _1,\mu _1)=\varphi (\beta _0,\mu _0)+\frac{1}{\beta }_{(\beta _0,\mu _0)}^{(\beta _1,\mu _1)}(E\mu x,\beta x)d(\beta ,\mu ),$$
(3)
where the integral is performed along a continuous path joining points $`(\beta _0,\mu _0)`$ and $`(\beta _1,\mu _1)`$ that does not encounter a phase transition. Semi-grand-canonical MC simulations can be directly used to determine the values of $`E`$ and $`x`$ for any given $`\beta `$ and $`\mu `$.

The potential at the initial point $`\varphi (\beta _0,\mu _0)`$ is usually chosen to be a convenient point where $`\varphi `$ can be computed analytically. One possible starting point is the high temperature limit ,

$$\beta \varphi (\beta ,\mu )\mathrm{ln}\left(2\right)\text{ as }\beta 0$$
(4)
where the limit is taken while keeping $`\mu `$ finite (implying that concentration $`x`$ converges to $`1/2`$ as $`\beta 0`$). Another possible starting point is in the limit of low temperature at a chemical potential stabilizing a given ground state $`g`$. In this limit, a “single-spin flip” low temperature expansion (LTE)  can be used to obtain:

$$\varphi (\beta ,\mu )=E_g\mu x_g\frac{1}{\beta N}\underset{s}{}\mathrm{exp}\left(\beta \left(\mathrm{\Delta }\epsilon _{s,g}\mu \mathrm{\Delta }\eta _{s,g}\right)\right),$$
(5)
where $`E_g`$ is the energy (per atom) of ground state $`g`$ with composition $`x`$, while $`\mathrm{\Delta }\epsilon _{s,g}`$ is the variation in the system’s total energy $`\left(NE\right)`$ associated with changing the identity of the atom sitting at site $`s`$ in ground state $`g`$, and $`\mathrm{\Delta }\eta _{s,g}`$ is the variation in $`\left(Nx\right)`$ associated with the same change. This infinite sum can be reduced to a finite number of terms by making use of the translational periodicity of the ground state.

In order to avoid lengthy MC simulations, it is useful to possess a criterion to find the highest possible temperature where the LTE is accurate. One attractive possibility is to compute $`\varphi `$ both from the LTE and from a mean field approximation (MFA). As long as both approaches give the same result within a user specified tolerance, we can be confident that “multiple spin flip events” are rare and that the LTE has the desired level of precision.

Using the LTE for each ground state of the system as a starting point to evaluate Equation (3), one can map out the whole potential surface $`\varphi ^\alpha (\beta ,\mu )`$ associated with any given ordered phase $`\alpha `$. Similarly, the high temperature limit can be used as a starting point to obtain $`\varphi ^\alpha (\beta ,\mu )`$ for the disordered phase. This process is illustrated in Figure 1. Once $`\varphi ^\alpha (\beta ,\mu )`$ has been determined for all phases, the boundary between two given phases $`\alpha `$ and $`\gamma `$ can be located by identifying the locus $`\mathrm{\Lambda }^{\alpha \gamma }=\{(\beta ,\mu ):\varphi ^\alpha (\beta ,\mu )=\varphi ^\gamma (\beta ,\mu )\}`$ where the potential surfaces intersect. At such intersections $`(\beta ,\mu )\mathrm{\Lambda }^{\alpha \gamma }`$, the concentrations $`x^\alpha `$ and $`x^\gamma `$ of the two phases in equilibrium are simply given by the slopes of the potential surface at the point of intersection:

$`x^\alpha `$ $`=`$ $`{\displaystyle \frac{\varphi ^\alpha (\beta ,\mu )}{\mu }}`$ (6)
$`x^\gamma `$ $`=`$ $`{\displaystyle \frac{\varphi ^\gamma (\beta ,\mu )}{\mu }}.`$ (7)
The calculations of thermodynamic quantities based on semi-grand-canonical MC simulations offer several advantages over their canonical counterparts.

1. For a given value of the control variables $`(\beta ,\mu )`$, the thermodynamic equilibrium of the system is never a phase-separated mixture, implying that the calculated quantities always reflect the properties of pure phases, free of interfacial contributions that would otherwise not be negligible due to the finite size of the simulation cell. (Note that the MC simulation cell must be commensurate with the periodicity of the equilibrium structure for this property to hold.)
2. The potential $`\varphi `$ can be very directly determined through a simple thermodynamic integration procedure (Equation (3)) where each required quantity takes the form of a thermodynamic average of computationally inexpensive quantities ($`E`$, which must already be computed in order to perform the simulation, and $`x`$, which is trivial to obtain).
3. Phase boundaries can be located by looking for intersections between curves, a criterion which is somewhat simpler to implement than the common tangent construction (although both methods are formally equivalent).

## 3 Algorithms

While the methods described in the previous section in principle enable the determination of thermodynamic quantities from MC simulations, a number of practical issues need to be addressed before this process can be fully automated.

A typical MC code requires a substantial amount of user intervention and relies on the user’s physical intuition to identify when a simulation is successful and reliable. The type of user interventions needed include:

1. Providing a variety of input parameters controlling the computational accuracy that need to be determined by trial and error.
2. Requiring the user to monitor the simulation periodically to ensure that the system has not undergone an unexpected phase transformation.
3. Performing a substantial amount of post-processing in order to obtain a phase diagram.

In this section we introduce algorithms to perform these tasks automatically and provide a formal justification of their applicability.

### 3.1 Reaching a target precision

Three parameters control the precision of a thermodynamic quantity computed through MC simulations: the simulation cell size, the equilibration time and the averaging time. While the analysis of the effect of simulation cell size on the precision is beyond the scope of the present paper, this section provides a sound basis for the selection of the two remaining parameters. The equilibration time is the time the system takes to reach thermodynamic equilibrium from a given initial configuration. Once equilibrium is reached, the instantaneous value of a given quantity (e.g. the energy) then needs to be averaged over a certain number of MC steps in order for its average to be sufficiently accurate. This defines the averaging time, which is the first parameter we will focus on.

#### 3.1.1 Averaging time

While the approach presented in this section relies on standard statistical results, it is nevertheless included for the purpose of introducing our notation and clarifying the assumptions made. Consider a microscopic quantity taking the value $`Q_t`$ at step $`t`$ and whose expectation value is equal to the desired thermodynamic quantity $`Q`$. If $`L`$ observations of $`Q_t`$ are available, a natural estimator of $`Q`$ is the average:

$$\overline{Q}_{[1,L]}=\frac{1}{L}\underset{t=1}{\overset{L}{}}Q_t.$$
(8)
The precision of this statistical quantity can be quantified by calculating its variance. In general, the variance of $`\overline{Q}_{[1,L]}`$ is given by:

$$Var\left(\overline{Q}_{[1,L]}\right)=\frac{1}{L^2}\underset{t=1}{\overset{L}{}}\underset{s=1}{\overset{L}{}}V_{\left|ts\right|}$$
(9)
where $`V_{\left|ts\right|}`$ is the covariance between $`Q_t`$ and $`Q_s`$. The fact that the sequence $`Q_t`$ is stationary when the system has reached thermodynamic equilibrium enables a simple estimator $`\widehat{V}_l`$ of the covariance structure $`V_l`$ of the sequence through the equation

$$\widehat{V}_l=\frac{1}{Ll}\underset{t=1}{\overset{Ll}{}}Q_tQ_{t+l}\overline{Q}_{[1,L]}^2.$$
(10)
While Equation 9 is very general, it is also very computationally intensive, requiring on the order of $`L^2`$ operations to estimate the variance of $`\overline{Q}_{[1,L]}`$. Following most of the literature on MC simulation, we consider the more computationally tractable special case obtained when the covariance structure $`V_l`$ is assumed to be of the form

$$V_l=V\rho ^{\left|l\right|}.$$
(11)
This approximation can be motivated as follows. A Markov chain can be represented by a linear superposition of components with exponentially decaying autocorrelation function. Each component is characterized by a certain correlation length and the component with the longest correlation length will typically give by far the largest contribution to the variance. Considering only the most slowly decaying component should thus provide a good approximation to the true variance of the quantity of interest.

Under the assumption that the averaging time is much larger than the longest correlation time,111This assumption should not be a concern because (i) this assumption tends to be violated only for averaging times such that the desired target precision on $`\overline{Q}`$ has not yet been reached, and (ii) our expression overestimates the true variance when this assumption is violated. “boundary” effects can be neglected and Equation (9) then reduces to

$$Var\left(\overline{Q}_{[1,L]}\right)=\frac{1}{L^2}\underset{t=1}{\overset{L}{}}\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}V\rho ^{\left|l\right|}.$$
(12)
Evaluating the geometric sum then yields:

$$Var\left(\overline{Q}_{[1,L]}\right)=\frac{V}{L}\left(\frac{1+\rho }{1\rho }\right).$$
(13)
Note that, since $`\rho =\mathrm{exp}\left(1/\tau \right)`$ where $`\tau `$ is the characteristic relaxation time, Equation (13) can be written, in the limit of large $`\tau `$, in the more familiar form $`Var\left(\overline{Q}_{[1,L]}\right)=2\tau V/L`$.

While $`V`$ can be directly estimated from Equation (10), setting $`l=0`$, the parameter $`\rho `$ can be found in a variety of ways. We found the following algorithm to be both fast and accurate: Find the time interval $`l^{}`$ such that $`\widehat{V}_l^{}/\widehat{V}_0=\frac{1}{2}`$ and then set $`\rho =2^{1/l^{}}`$. This approach offers the important advantage that it uses, as an input, only the correlation between relatively distant values of $`Q_t`$. This ensures that $`\rho `$ will be determined by the most slowly decaying components of the Markov process, which are the ones that the functional form $`V\rho ^{\left|l\right|}`$ was meant to model. It is also a very fast algorithm, since a simple bracketing algorithm allows the determination of $`\rho `$ in about $`L\mathrm{log}L`$ operations. In the event that there are more than one $`l^{}`$ satisfying the equality, we chose the smallest one. The rationale behind this choice is that there is a small but nonzero chance that the estimated covariances are such that $`\widehat{V}_l/\widehat{V}_0=1/2`$ although the true covariances are *not* such that $`V_l/V=1/2`$. Choosing the smallest $`l^{}`$ makes the procedure less sensitive to these spurious events since the probability of occurence of one such event increases with the range of values of $`l`$ considered.

In short, the averaging time $`L`$ needed to reach a given precision can be found by computing the variance of $`\overline{Q}_{[1,L]}`$ through Equation (13) periodically and stopping the simulation when the variance goes below a user specified threshold:

$$Var\left(\overline{Q}_{[1,L]}\right)\left(\frac{p}{z_\alpha }\right)^2,$$
(14)
where $`p`$ is the target precision on the value of $`Q`$ while

$$z_\alpha =\sqrt{2}\text{erf}^1\left(1\alpha \right)$$
(15)
is the critical value associated with the desired level of confidence $`\alpha `$. For instance, if the error on $`Q`$ is required to be less than $`10^3`$ with a probability of $`99\%`$, then $`p=10^3`$, $`\alpha =0.01`$ and $`z_\alpha =2.576`$.

#### 3.1.2 Equilibration time

When devising a criterion to determine whether a simulation has reached thermodynamic equilibrium, approaches based on identifying correlation lengths can be unreliable because the system may be very far from equilibrium during equilibration. For this reason, we employ a relatively simple, and yet very robust approach that is based on the main defining characteristic of thermodynamic equilibrium. If the average of the quantity $`Q_t`$ between steps $`t_1`$ and $`t_2`$ is not statistically significantly different from the average of $`Q_t`$ between steps $`t_2`$ and $`t_3`$, then the hypothesis that equilibrium has been reached at step $`t_1`$ and beyond cannot be rejected. Of course, the power of this test increases as $`t_2t_1`$ and $`t_3t_2`$ increase and this definition of equilibrium thus depends on the target precision $`p`$ we seek to reach. This is to be expected since, strictly speaking, thermodynamic equilibrium is never reached in a simulation of a finite length. Any definition of equilibrium must include an arbitrary cutoff and it is natural to choose this cutoff to be a function of the precision we require.

Both the equilibration and averaging time needed to reach a precision $`p`$ can thus be determined through an algorithm that can be outlined as follows (refer to Figure 2). Let $`L`$ denote the number of samples of $`Q_t`$ collected so far. Consider a range of values of $`t_1[1,L]`$ and for each, check whether

$$\mathrm{\Delta }\overline{Q}\left|\overline{Q}_{[t_1,\left(t_1+L\right)/2]}\overline{Q}_{[\left(t_1+L\right)/2,L]}\right|p,$$
(16)
where $`Q_{[t_1,t_2]}`$ denotes the average of $`Q_t`$ for $`t[t_1,t_2]`$. Then, verify that $`Var\left(\overline{Q}_{[t_1,L]}\right)`$ satisfies Equation (14). If both conditions are satisfied, then $`\overline{Q}_{[t_1,L]}`$ provides an estimate of $`Q`$ with a precision $`p`$ and with a confidence level $`\alpha `$.

In its simplest form, the implementation of this algorithm involves storing *all* samples of $`Q_t`$ collected so far and requires a loop over values of $`t_1`$ which make the complexity of the algorithm $`O\left(L^2\right)`$. Various techniques help reduce the resources required. First, sampling $`Q_t`$ every MC pass,222A MC pass consists of a number of attempted spin flips equal to the number of lattice sites in the simulation cell. instead of at every spin flip, results in essentially no loss in precision since the correlation length of the simulation typically extends much beyond a MC pass. Second, one can limit the search for $`t_1`$ to values that are integer multiple of some block size $`L_b`$ that is allowed to increase linearly with the total number of MC passes $`L`$. In this fashion, the number of operations needed to determine $`t_1`$ does not increase as the simulation progresses. The value of $`t_1`$ thus found by a search constrained to block boundaries differs at most by $`L_b`$ from the $`t_1`$ which would have been found without block constraints. In the worst case, the simulation may have to be run for an additional $`L_b`$ passes and it follows that the length of the simulation will be at most $`\left(1+L_b/L\right)`$ larger than the optimal computational time that would have been needed without block constraints. A computationally convenient way to implement this scheme is to multiply the block size $`L_b`$ by two whenever the total number of passes $`L`$ doubles. Thanks to this approach, the entire history of values of $`Q_t`$ does not need to be stored, only block averages of $`\overline{Q}`$, $`V`$ and $`\rho `$ are needed, as illustrated in Figure 3. When the block size doubles, the values of $`\overline{Q}`$, $`V`$ and $`\rho `$ for two consecutive blocks can be combined to obtain the corresponding value for the new block having twice the length.333The only drawback is that the resulting estimate of $`\rho `$ does not make use of all the available data points, because the product $`Q_tQ_s`$ for $`t`$ and $`s`$ belonging to different blocks cannot be determined from within-block averages. Although it is inefficient, our estimate of $`\rho `$ is nevertheless unbiased and the loss of efficiency is negligible if the block size is large relative to the correlation time.

### 3.2 Detecting phase transitions

Phase transitions can be identified by locating singularities in some function $`Q\left(C\right)`$, where $`Q`$ is some thermodynamic quantity (such as energy, concentration or an order parameter) and $`C`$ is some control variable (such as temperature or chemical potential). In this section, we propose a simple algorithm to perform the detection of such singularities while accounting for the fact that quantities obtained from MC are intrinsically noisy. Our algorithm is applicable to the detection of both first-order and second-order transitions, although the power of the test is clearly higher for first-order transitions.

An analytic function has the property that the knowledge of its shape in some finite interval enables the prediction of its value outside of that interval. The failure of this extrapolation procedure in the presence of a singularity permits its detection. Formally, an analytic function can be extrapolated using a Taylor series. We approximate this operation by fitting the output of the MC simulation by a polynomial.

Consider $`n+1`$ estimates $`\overline{Q}_1,\mathrm{},\overline{Q}_{n+1}`$ of the quantity $`Q\left(C\right)`$ obtained for a sequence of values of the control variable $`C_1,\mathrm{},C_{n+1}`$. The general principle behind our singularity detection procedure is to fit a polynomial through the first $`n`$ points and if the observed value $`\overline{Q}_{n+1}`$ is far from the value predicted from the polynomial fitted to the previous data, a singularity must be present between point $`n`$ and $`n+1`$ (see Figure 4). To make this algorithm practical, two questions need to be answered. First, how do we choose the order of the polynomial to be fitted? Second, how do we select the maximum allowable prediction error that can be tolerated without claiming the presence of a singularity?

To answer the first question, we make use of the well-known cross-validation criterion : Choose the order $`k`$ of the polynomial that minimizes the quantity:

$$CV=\frac{1}{n}\underset{i=1}{\overset{n}{}}\left(\overline{Q}_iQ_i^{}\right)^2$$
(17)
where $`Q_i^{}`$ is the value of $`\overline{Q}_i`$ predicted from a polynomial least-squares fit to the remaining $`n1`$ points: $`\overline{Q}_1,\mathrm{},\overline{Q}_{i1},\overline{Q}_{i+1},\mathrm{},\overline{Q}_n`$. Intuitively, this criterion tests the predictive power of the polynomial using the points $`\overline{Q}_1,\mathrm{},\overline{Q}_n`$ before we attempt to use it for the prediction of $`\overline{Q}_{n+1}`$. Once the order $`k`$ of the polynomial has been chosen, the standard theory of least-squares  provides the distribution of the predicted value $`Q_{n+1}^{}`$, which can be used to determine whether the prediction error is statistically significantly different from zero. For $`i=1,\mathrm{},n`$ and $`j=1,\mathrm{},k+1`$, let

$`X_{ij}`$ $`=`$ $`\left(C_i\right)^{j1}`$ (18)
$`Y_i`$ $`=`$ $`\overline{Q}_i,`$ (19)
so that the equation defining the least-squares problem can be written in matrix form as $`Y=Xa+\epsilon `$, where $`\epsilon `$ is an $`n\times 1`$ vector of error terms. The vector $`a`$ of the coefficients of the polynomial can be estimated by

$$\widehat{a}=\left(X^TX\right)^1X^TY$$
(20)
while the covariance matrix of $`\widehat{a}`$ is given by

$$V=\sigma ^2\left(X^TX\right)^1$$
(21)
where $`\sigma ^2`$ is the variance of the residuals $`\epsilon _i`$ of the fit. The quantity $`\sigma ^2`$ can be estimated in two asymptotically equivalent ways. Equation (13) provides an estimate of the variance of any of the $`\overline{Q}_i`$, which can be averaged over $`i`$ to yield a single, more accurate estimate:

$$\widehat{\sigma }^2=\frac{1}{n}\underset{i=1}{\overset{n}{}}Var\left(\overline{Q}_i\right).$$
(22)
Alternatively, the residuals of the least-squares fit can be used:

$$\widehat{\sigma }^2=YXa^2/\left(nk1\right).$$
(23)
While Equation (22) is more accurate (especially if $`n`$ is small), because it makes use of additional information, Equation (23) is simpler to implement as it minimizes the amount of information that needs to be transferred between the Monte Carlo simulation code and the phase transition detection routine.

The predicted value of $`Q`$ for point $`n+1`$ is then given by

$$Q_{n+1}^{}=x^T\widehat{a}$$
(24)
where $`x_j=\left(C_{n+1}\right)^{j1}`$ for $`j=1,\mathrm{},k+1`$, while the variance of this prediction is given by:

$$v=x^TVx.$$
(25)
If there is no singularity, the variance of the difference between $`\overline{Q}_{n+1}`$ and $`Q_{n+1}^{}`$ is $`v+\sigma ^2`$because (i) the variance of $`Q_{n+1}^{}`$ is $`v`$ while the variance of $`\overline{Q}_{n+1}`$ is $`\sigma ^2`$, (ii) $`\overline{Q}_{n+1}`$ and $`Q_{n+1}^{}`$ are independent. An asymptotically valid statistical test thus consists in rejecting the hypothesis that there is no phase transition if the prediction error $`\left|\overline{Q}_{n+1}Q_{n+1}^{}\right|`$ is such that444Note that $`Q_{n+1}^{}`$ must be an unbiased predictor of $`Q\left(C_{n+1}\right)`$ for this test to be valid, which is not the case for finite $`n`$ due to the finite order of the approximating polynomial. Fortunately, the bias (squared) goes to zero as $`n\mathrm{}`$ and eventually becomes negligible relative to $`\sigma ^2`$, ensuring the asymptotic validity of the procedure. Note that the variance $`v`$ of the prediction is also negligible asymptotically relative to $`\sigma ^2`$. We nevertheless include the variance term $`v`$ in Equation (26) in order to improve the accuracy of the test for finite $`n`$, since we do have an estimator of the variance (unlike the bias).

$$\left|\overline{Q}_{n+1}Q_{n+1}^{}\right|z_\alpha \sqrt{v+\sigma ^2}.$$
(26)
For $`z_\alpha `$ defined as in Equation (15), the probability that a phase transition is incorrectly identified is less than $`\alpha `$. The relevant quantities are represented in Figure 4.

The Equations presented above rely on the assumption that the deviations away from the true thermodynamic quantity, $`\left(\overline{Q}_iQ\left(C_i\right)\right)`$ for $`i=1,\mathrm{},n+1`$, are statistically independent. This is an appropriate assumption if the system is given sufficient time to equilibrate every time the value of the control variable $`C`$ changes. Our method also assumes that the deviations $`\left(\overline{Q}_iQ\left(C_i\right)\right)`$ all have the same variance. Both of these requirements are automatically satisfied if the points $`\overline{Q}_1,\mathrm{},\overline{Q}_{n+1}`$ are generated using the algorithm presented in Section 3.1.

For the sole purpose of locating phase transitions, one particular type of thermodynamic function enables an especially sensitive test for phase transitions: order parameters. Unfortunately, order parameters that are invariant under any possible symmetry operation of the lattice are often expensive to compute. However, if one always traverses phase transitions from an ordered phase toward a disordered phase (or another ordered phase), the task is dramatically simplified. When a MC cell is initialized with one specific ground state (i.e. an ordered phase) and before a phase transition is reached, there is a vanishing probability that the system will fluctuate enough to reach a state that is different but symmetrically equivalent to the original state. Thus, an order parameter can be computed by simply keeping track of the number of sites that host an atomic specie that is different from the specie found at the same site at initialization. There is no need to check for every possible symmetry-related ordered configuration.

The algorithms just presented are especially useful when one desires to construct the potential surface $`\varphi (\beta ,\mu )`$ of a phase (using thermodynamic integration over the paths depicted in Figure 1) without requiring the user to manually determine when to stop the thermodynamic integration along a given direction because a phase transition has occurred. The automated construction of such free energy surfaces provides a very natural pathway through which an atomistic MC simulation can be used to provide input to CALPHAD-type calculations based on phenomenological free energy models.

### 3.3 Phase boundary tracing

When one is solely interested in determining a system’s phase diagram, it would be computationally advantageous to be able to follow the boundary of a phase without first calculating the potential $`\varphi `$ over the whole region where this phase is stable. In this section, we describe how, in the case of a first-order transition, the entire boundary of a two-phase equilibrium can be determined, starting from a single point where a transition is known to occur. We first focus on a simple method that neglects the unavoidable presence of statistical noise in MC simulations before describing how the presence of this noise can be handled.

Consider a value of the reciprocal temperature $`\beta `$ and of the chemical potential $`\mu `$ that is known to stabilize a phase-separated mixture between phases $`\alpha `$ and $`\gamma `$. The knowledge of such a pair $`(\beta ,\mu )`$ could be provided, for instance, by the phase transition detection algorithm presented in the previous section, assuming that the hysteresis loop of the phase transition is sufficiently narrow (which tends to be case at high temperature). Alternatively, if $`T=0`$, the chemical potential that simultaneously stabilizes two ground states $`\alpha `$ and $`\gamma `$ can be found by a simple common tangent construction from the knowledge of the energy and composition of each ground state. At such a two-phase equilibrium point, the potentials $`\varphi `$ of each phase are equal. Our task is to find, as reciprocal temperature $`\beta `$ changes by an amount $`d\beta `$, the changes in chemical potential needed to preserve the equality between each phase’s thermodynamic potential. By iterating this process, one can define a whole path $`\mu \left(\beta \right)`$ that traces out the value of the chemical potential at the phase transition as temperature changes, and thus find the concentrations $`x_\alpha `$ and $`x_\gamma `$ of each phase in equilibrium as a function of temperature. The differential equation whose solution provides the desired path $`\mu \left(\beta \right)`$ can be found by taking the total differential of the equation $`\beta \varphi ^\alpha =\beta \varphi ^\gamma `$ following from the identity between the potentials $`\varphi `$ of each phase:

$`{\displaystyle \frac{\left(\beta \varphi ^\alpha \right)}{\beta }}|_\mu d\beta +{\displaystyle \frac{\left(\beta \varphi ^\alpha \right)}{\mu }}|_\beta d\mu ={\displaystyle \frac{\left(\beta \varphi ^\gamma \right)}{\beta }}|_\mu d\beta +{\displaystyle \frac{\left(\beta \varphi ^\gamma \right)}{\mu }}|_\beta d\mu `$ (27)
$`\left(E^\alpha \mu x^\alpha \right)d\beta \beta x^\alpha d\mu =\left(E^\gamma \mu x^\gamma \right)d\beta \beta x^\gamma d\mu `$ (28)
$`{\displaystyle \frac{d\mu }{d\beta }}={\displaystyle \frac{\left(E^\gamma E^\alpha \right)}{\beta \left(x^\gamma x^\alpha \right)}}{\displaystyle \frac{\mu }{\beta }},`$ (29)
where $`E^\alpha `$ and $`x^\alpha `$ are the energy and concentration of phase $`\alpha `$ at reciprocal temperature $`\beta `$ and chemical potential $`\mu `$ (and similarly for phase $`\gamma `$). All the quantities needed can be directly obtained from MC simulations. Thanks to the fact that both phases are metastable for values of the chemical potential in the neighborhood of the true phase transition, it is possible to obtain $`E`$ and $`x`$ for both phases at the same chemical potential. A very efficient way to implement this algorithm is to keep two simulations in memory at the same time — one for the phase $`\alpha `$ and one for the phase $`\gamma `$. In this fashion, the final configuration of the simulation of phase $`\alpha `$ at reciprocal temperature $`\beta `$ can be used as an initial configuration for the simulation of phase $`\alpha `$ at temperature $`\beta +d\beta `$ (and similarly for phase $`\gamma `$). The long equilibration time associated with the crossing of a phase transition is thus avoided. Interestingly, this algorithm can be directly used even when the two phases in equilibrium are based on a completely different parent lattice.

We now address the issue that the quantities obtained from MC are contaminated by noise. The presence of noise causes the solution to Equation (29) calculated from MC to deviate from the true path $`\mu \left(\beta \right)`$. This eventuality brings up two questions. First, is the solution to Equation (29) stable or unstable with respect to the presence of small perturbations? Second, what can one do if the noise becomes so large that the calculated chemical potential leaves the region where both phases are metastable?

To answer the first question, consider a perturbation $`\epsilon \left(\beta \right)`$ of the chemical potential away from its value $`\mu \left(\beta \right)`$ at the phase transition at reciprocal temperature $`\beta `$. We now calculate how this error is further propagated when solving Equation 29. To the first order, the difference between the potential $`\varphi `$ of the two phases is $`\beta \left(x^\gamma \left(\beta \right)x^\alpha \left(\beta \right)\right)\epsilon \left(\beta \right)`$. As temperature changes, Equation (29) updates $`\mu `$ so as to preserve this difference. We thus have the equality:

$$\left(x^\gamma \left(\beta \right)x^\alpha \left(\beta \right)\right)\epsilon \left(\beta \right)=\left(x^\gamma \left(\beta +d\beta \right)x^\alpha \left(\beta +d\beta \right)\right)\epsilon \left(\beta +d\beta \right).$$
(30)
Taking the limit of infinitesimal $`d\beta `$ and rearranging yields:

$$\frac{d\epsilon }{d\beta }=\epsilon \frac{d}{d\beta }\mathrm{ln}\left|x^\gamma x^\alpha \right|$$
(31)
and it follows that the small perturbation $`\epsilon `$ is unstable as $`\beta `$ increases when the difference in concentration between the two phases in equilibrium $`\left|x^\gamma x^\alpha \right|`$ decreases with increasing $`\beta `$. The practical consequence of this observation is that integrating Equation (29) from high temperature to low temperature yields a method more robust to statistical noise, because $`\left|x^\gamma x^\alpha \right|`$ usually increases with decreasing temperature $`T`$.

We now address the second robustness issue, namely, what can we do if the algorithm wanders outside of the region of metastability of either phase? It is straightforward to detect when this happens by simply checking whether $`\left|x^\gamma x^\alpha \right|`$ is not statistically significantly different from zero. One then needs to determine which of the two MC simulations (for phase $`\alpha `$ or $`\gamma `$) underwent a phase transition. This can be accomplished by keeping track of the concentration of each phase ($`x_1^\alpha ,\mathrm{},x_n^\alpha `$ and $`x_1^\gamma ,\mathrm{},x_n^\gamma `$) at temperatures previously visited ($`\beta _1,\mathrm{},\beta _n`$) and using them the verify whether the new concentration $`x_{n+1}^\alpha x_{n+1}^\gamma `$ at $`\beta _{n+1}`$ is best predicted by the sequence $`x_i^\alpha `$ or $`x_i^\gamma `$, using the extrapolation described in the previous section (see Figure 5). The sequence that gives the worst prediction will point to the MC simulation that underwent a phase transition.

Without loss of generality consider the case when $`x^\alpha <x^\gamma `$ and when the simulation of phase $`\alpha `$ transformed to phase $`\gamma `$. We can then “recenter” the chemical potential to a good estimate of its true value at the phase transition as follows (refer to Figure 6).

1. Gradually decrease the chemical potential of the MC simulation of phase $`\alpha `$ until its concentration can again be predicted by the sequence $`x_i^\alpha `$ within the magnitude of the statistical noise. Let $`\underset{¯}{\mu }`$ denote the value of the chemical potential when this happens and save the state $`S`$ of simulation cell at that point.
2. Then, gradually increase the chemical potential of the MC simulation of phase $`\alpha `$ until its concentration can be predicted by the sequence $`x_i^\gamma `$ within the magnitude of the statistical noise. Let $`\overline{\mu }`$ denote the value of the chemical potential when this happens.
3. Restore the state of the simulation cell of phase $`\alpha `$ to the saved state $`S`$ and set the chemical potential to $`\mu =\left(\overline{\mu }+\underset{¯}{\mu }\right)/2`$.
4. Rerun the simulation for both phases $`\alpha `$ and $`\gamma `$ at chemical potential $`\mu `$ and continue the integration Equation (29) with the values of $`E`$ and $`x`$ obtained for each phase.

This approach results in a very robust algorithm that can run without user intervention to map out the whole two-phase equilibrium. Thanks to this safety mechanism, we found that it is feasible to integrate Equation (29) towards increasing temperature, despite the potential for instability with respect to noise. This proves useful since it is very easy to find a starting point for our algorithm at $`T=0`$. The boundary tracing algorithm may have to recenter itself a few times on its way to the end of the two-phase equilibrium. But once the location of the phase boundary at high temperature is known, a smooth phase boundary can then be calculated, starting from the highest temperature and following the phase boundary in the direction of decreasing temperature.

The end of the two-phase region can be detected by monitoring the occurrence of either one the following events (see Figure 7):

1. The concentrations $`x^\alpha `$ and $`x^\gamma `$ are statistically significantly different, but either $`x^\alpha `$ or $`x^\gamma `$ is not predicted by the corresponding sequence $`x_i^\alpha `$ or $`x_i^\gamma `$within statistical uncertainty.
2. The concentrations $`x^\alpha `$ and $`x^\gamma `$ are not statistically significantly different, but the sequences $`x_i^\alpha `$ and $`x_i^\gamma `$ both predict the same concentration, within statistical uncertainty.

Case 1 indicates that a third phase $`\delta `$ has become stable (i.e. a eutectic or peritectic point), while Case 2 indicates that the phase boundaries have met at a common composition (e.g. at a congruent point or a critical point). In the first case, the algorithm will actually continue slightly beyond the true eutectic or peritectic point because of the presence of hysteresis at a first-order phase transition. However, the true location of the three-phase equilibrium can be determined from the intersection of the $`\alpha \delta `$ and $`\delta \gamma `$ phase boundaries, once they have been determined.

## 4 Applications

This section presents a number of examples of automated calculations of thermodynamic properties using the previously introduced algorithms. In each example, the cluster expansion describing the energetics of the alloy system was obtained from a fit to first-principles structural energy calculations using the MAPS  code, which automates the process. A description of the algorithms enabling the automatic construction of such a cluster expansion, as well as a description of the interaction parameters defining each cluster expansion used below, can be found in .

Our first example is the determination of the complete potential surface $`\varphi (T,\mu )`$ of the solid-state phases of the MgO-CaO pseudo-binary alloy system (see Figure 8). As the MgO-CaO system is phase separating (as shown in Figure 12), the thermodynamic potential surface exhibits a characteristic butterfly shape. Below the critical temperature, the surface intersects with itself at a nonzero angle, indicating the presence of a miscibility gap. Since the MC simulation is able to sample metastable phases, for some values of the chemical potential, one can determine $`\varphi (T,\mu )`$ for both phases in the region of metastability, which explains the double-valued nature of the potential surface at low temperatures. Above the critical temperature, $`\varphi (T,\mu )`$ becomes single-valued and varies smoothly as $`\mu `$ varies. (Note that, in reality, the alloy melts before reaching that point at atmospheric pressure.) Figure 9 shows the Gibbs555The difference between the Gibbs and Helmholtz free energies can be neglected at atmospheric pressure in solid state systems because the product of pressure $`P`$ with changes in specific volume ($`\mathrm{\Delta }V`$) are small. free energy surface $`G(x,T)\varphi (T,\mu (T,x))+\mu (T,x)x`$ for the same system and, in this more familiar representation, the presence of miscibility gap at low temperature is clearly recognizable. The calculated Gibbs free energy surface can be directly used to fit standard solution models, such as an ideal solution model, plus a polynomial in temperature and concentration to account for non-ideality. The order of the polynomial needed to represent the MC output can, for instance, be determined through cross-validation, as defined in Section 3.2. The resulting very compact representation of the thermodynamics of the alloy system can be used in CALPHAD-type calculations to supplement existing experimental data or to fill in data not yet determined experimentally.

Figure 10 shows the result of a similar calculation in the case of a ordering alloy system, the Ti-Al system. For clarity, the metastable portions of the $`\varphi (T,\mu )`$ surface are omitted so that the potential is single-valued everywhere. The order-disorder transition is clearly visible from the kink in the surface, which is highlighted by a thick curve. The corresponding Gibbs free energy representation of the same data is shown in Figure 11.

Figure 12 illustrates the usefulness of the boundary tracing algorithm described in Section 3.3. The examples shown even include a case (Ti-Al) where a equilibrium between two phases differing by their lattice type: the L10 is fcc-based while the DO19 and the Ti solid solution are hcp-based). These examples show that one can obtain phase diagrams from first-principles without first obtaining the complete $`\varphi (T,\mu )`$ surface.

It is important to note that, in these examples, the MC code requires very little intervention or input from the user. Thanks to the algorithms introduced in Section 3.1, the code is able to autonomously determine, at each value of the control variables $`(T,\mu )`$, when equilibrium has been reached and how long the simulation needs to be run in order to reach the target precision. The code is also able to determine when the region of stability of one phase has been exited, thanks to the algorithms described in Section 3.2, so that the user does not need to specify when to stop the thermodynamic integration. The only input parameters the user needs to specify are

1. which phase to scan (the code can determine the value of the chemical potential $`\mu `$ stabilizing the desired two-phase equilibrium at $`0`$K);
2. the target precision (here, $`0.1`$ at % precision on the concentration variable is used as a criterion);
3. the size of the simulation cell (in the present study, the simulation cell is always chosen such that it contains a sphere of a diameter equal to at least 10 times the nearest-neighbor interatomic spacing).

## 5 Conclusion

The algorithms we have introduced enable researchers to focus on high-level concepts when employing MC simulations. All the tasks that make traditional Monte Carlo simulations so tedious have been formalized into algorithms that can autonomously initialize and control the appropriate MC simulation. The user only needs to specify high-level goals, such as the target precision or which phase to focus on. This automation of MC simulations, along with the earlier work that enabled the automation of cluster expansion construction , signifies that it is now possible for large community of researchers to employ first-principles calculations to augment and complete existing experimental thermodynamic databases without facing the steep learning curve that has traditionally been associated with first-principles thermodynamical calculations.

This work was supported by the NSF under program DMR-0080766.

## References


format: text

text:
# Electromigration-Induced Propagation of Nonlinear Surface Waves

## Abstract

Due to the effects of surface electromigration, waves can propagate over the free surface of a current-carrying metallic or semiconducting film of thickness $`h_0`$. In this paper, waves of finite amplitude, and slow modulations of these waves, are studied. Periodic wave trains of finite amplitude are found, as well as their dispersion relation. If the film material is isotropic, a wave train with wavelength $`\lambda `$ is unstable if $`\lambda /h_0<3.9027\mathrm{}`$, and is otherwise marginally stable. The equation of motion for slow modulations of a finite amplitude, periodic wave train is shown to be the nonlinear Schrödinger equation. As a result, envelope solitons can travel over the film’s surface.

I. Introduction

Waves propagating over the free surface of a layer of water subject to gravity have been studied extensively for well over a century. In his ground breaking paper of 1847, Stokes showed that steady periodic wave trains of finite amplitude are possible, and that their dispersion relation depends upon their amplitude Stokes . In these wave trains, the effects of dispersion and nonlinearity balance, yielding a steady wave form.

In the 1960’s, it was discovered that if the water is of finite depth $`h_0`$ and the wavelength $`\lambda =2\pi /k`$ is sufficiently short, the wave train is unstable. Benjamin and Feir BF  showed that the Stokes wave train is unstable to small modulational or side band perturbations if $`kh_0>1.363\mathrm{}`$, and this been confirmed experimentally Benjamin ; Yuen ; Lake ; Shemer .

Benjamin and Feir’s analysis holds only at early times, just after the instability has begun to develop. Zakharov Zakharov  and later Hasimoto and Ono HO  derived an equation of motion for the envelope of a packet of finite-amplitude gravity waves which remains valid long after the onset of the instability. Their equation — the nonlinear Schrödinger equation — admits envelope soliton solutions that are intimately related to the ultimate state of the fluid surface.

Although it is not yet widely appreciated, there are nontrivial electrical free boundary problems. When an electrical current passes through a piece of solid metal or semiconductor, collisions between the conduction electrons and the atoms at the surface lead to drift of these atoms. This phenomenon, which is known as surface electromigration (SEM), can cause the surface of a solid to move and deform Latyshev ; Stoyanov ; Krug\_and\_Dobbs ; Kraftapl ; Suojap ; Mohanjap ; Krug1 ; Krug2 ; Gengor1 ; Gengor2 ; Gengor3 ; Mohanepl ; Mohanpre . The free surface of a conducting film moves in response to the electrical current flowing through the bulk of the film, in much the same way that flow in the bulk of a fluid affects the motion of its surface. However, the analogy is not perfect — the boundary conditions are very different in the two problems. In addition, metals and semiconductors are anisotropic materials, and this distinguishes them from fluids as well.

Waves can propagate over the free surface of a current-carrying film deposited on an insulating substrate, at least in the limit of small amplitude. These waves, whose propagation is induced by surface electromigration, are somewhat analogous to gravity waves, and will be called “SEM waves.” The linear dispersion relation for these waves has been known for some time Krug1 , and will be discussed briefly in Sec. III.

It is natural to ask to what extent the nonlinear behavior of gravity waves is paralleled by that of SEM waves. To address this question, in this paper I will study SEM waves of finite amplitude, and slow modulations of these waves. Periodic wave trains of finite amplitude (the analog of the Stokes wave) are found, as well as their dispersion relation. If the film material is isotropic, the wave trains are unstable for $`kh_0>1.609927\mathrm{}`$ and are otherwise marginally stable. The equation of motion for slow modulations of the finite amplitude, periodic wave train is shown to be the nonlinear Schrödinger equation.

It is worth mentioning that SEM is not just of academic interest: In modern integrated circuits, the current-carrying metal lines or “interconnects” carry very high current densities, and electrical failure can occur due to the effects of SEM. In particular, SEM can cause a small perturbation at the edge of a metal strip to become a slit-shaped void which propagates across the line until its tip contacts the opposite side of the strip; electrical failure is then complete Mohanepl ; Mohanpre ; Rose ; Sanchezapl .

II. Equations of Motion

Consider a single-crystal metallic or semiconducting film of uniform thickness $`h_0`$ deposited on the plane surface of an insulating substrate. We take the $`z`$ axis to be normal to the film-vacuum interface and locate the origin $`O`$ in this plane. (We will occasionally find it convenient to use another set of Cartesian coordinates $`(x^{},z^{})`$ with $`x^{}=x`$, $`z^{}=z+h_0`$ and origin $`O^{}`$ at the substrate surface.) The film’s surface $`z=0`$ will be assumed to be a low index crystal plane. A constant electrical current flows through the film in the $`x`$ direction, and the electric field within the film is $`\stackrel{}{E}_0=E_0\widehat{x}`$.

Now suppose that the upper surface of the film is perturbed. Let the outward-pointing unit normal to this surface be $`\widehat{n}`$. For simplicity, we shall restrict our attention to perturbations whose form does not depend upon $`y`$, so that the height of the film’s surface above the substrate $`h=h_0+\zeta `$ depends only on $`x`$ and $`t`$ (Fig. 1). The upper film surface will evolve in the course of time due to the effects of SEM and surface self-diffusion. We assume that the current flowing through the film is held fixed.

Clearly, the problem is two-dimensional (2D), and the dependence of all quantities on $`y`$ will therefore be suppressed. The electrical potential $`\mathrm{\Phi }=\mathrm{\Phi }(x,z,t)`$ satisfies the 2D Laplace equation

$$^2\mathrm{\Phi }=0$$
(1)
and is subject to the boundary condition $`\widehat{n}\stackrel{}{}\mathrm{\Phi }=0`$ on the upper surface and $`\widehat{z}\stackrel{}{}\mathrm{\Phi }=0`$ on the lower. More explicitly, we have

$$\mathrm{\Phi }_z(x,\zeta ,t)=\zeta _x\mathrm{\Phi }_x(x,\zeta ,t)$$
(2)
and

$$\mathrm{\Phi }_z(x,h_0,t)=0,$$
(3)
where $`f_xf/x`$ and so forth. If the initial perturbation is localized, we also have

$$\zeta 0$$
(4)
for $`x\pm \mathrm{}`$ and, furthermore,

$$\stackrel{}{}\mathrm{\Phi }(x,z,t)\stackrel{}{E}_0$$
(5)
for $`x\pm \mathrm{}`$ and $`h_0z0`$.

We assume that the atomic mobility is negligible at the film-insulator interface, so that the form of that interface remains planar for all time. Further, in the interest of simplicity, we assume that the applied current is high enough that the effects of SEM are much more important than those of capillarity. The surface atomic current $`\stackrel{}{J}_a`$ is then proportional to the electrical current at the surface. Explicitly, $`\stackrel{}{J}_a=\nu Mq\stackrel{}{}\mathrm{\Phi }`$, where the $`\nu `$ is the areal density of the mobile surface atoms, $`q`$ is their effective charge, and $`M`$, the adatom mobility, in general depends on the surface slope $`\zeta _x`$. When there is a net influx of atoms into a small surface element, it will move. The normal velocity of the film’s free surface is $`v_n=q\nu \mathrm{\Omega }_s(M_s\mathrm{\Phi })`$, where $`\mathrm{\Omega }`$ is the atomic volume and $`s`$ is the arc length along the surface. Therefore

$$\zeta _t=q\nu \mathrm{\Omega }_x\left[\frac{M(\zeta _x)}{\sqrt{1+\zeta _x^2}}_x\mathrm{\Phi }(x,\zeta ,t)\right],$$
(6)
where $`_x`$ is the total derivative with respect to $`x`$. (Since $`\zeta `$ depends on $`x`$, we have $`_x\mathrm{\Phi }(x,\zeta ,t)=\mathrm{\Phi }_x(x,\zeta ,t)+\zeta _x\mathrm{\Phi }_z(x,\zeta ,t)`$, for example.) Together, Eqs. (1) - (6) completely describe the nonlinear dynamics of the film surface.

III. Finite Amplitude Periodic Wave Trains

We wish to study the propagation of a disturbance whose amplitude is small. To this end, we put $`\mathrm{\Phi }=E_0x+\varphi `$ and presume that $`\zeta `$ and $`\varphi `$ are of order $`ϵ`$.

We shall suppose that the crystal structure is invariant under the reflection $`xx`$, so that $`M`$ is an even function of the slope $`\zeta _x`$. As a result, $`M(\zeta _x)=M_0(1+c\zeta _x^2)+O(ϵ^4)`$, where $`M_0M(0)`$ and $`c`$ is a dimensionless constant which vanishes for an isotropic material.

If we work to order $`ϵ`$, the equations of motion are linearized and we find that sinusoidal waves of the form $`\zeta =2a\mathrm{cos}(kx\omega _0t)`$ propagate over the surface. $`\omega _0=\omega _0(k)`$ is given by the linear dispersion relation

$$\omega _0=\beta E_0\frac{k^2}{\sigma },$$
(7)
where $`\beta qM_0\nu \mathrm{\Omega }`$ and $`\sigma \mathrm{tanh}kh_0`$ Krug1 . The corresponding group velocity is

$$v_g=v_0\left[h_0^2k^2\left(1\frac{1}{\sigma ^2}\right)+2\frac{h_0k}{\sigma }\right],$$
(8)
where $`v_0=\beta E_0/h_0`$ is the phase velocity of waves of infinitely long wavelength. If the amplitude $`2a`$ is not a constant but instead varies slowly with position, the resulting amplitude modulation will propagate with the group velocity.

The goal of this paper is to go beyond the linear approximation and study wave trains of finite amplitude. To begin, we will find the analog of the Stokes wave propagating over deep water: we will consider the limit $`h_0\mathrm{}`$ and restrict our attention to periodic wave trains with constant amplitude $`a`$. These simplifying assumptions will be relaxed in the next section.

Expanding $`\zeta `$ and $`\varphi `$ in powers of the amplitude, we have

$$\zeta =[a+\mu _2a^2^2+\mu _3a^3^3+\mathrm{}]+\mathrm{c}.\mathrm{c}.$$
(9)
and

$$\varphi =[\nu _1ae^{kz}+\nu _2a^2e^{2kz}^2+\nu _3a^3e^{3kz}^3+\stackrel{~}{\nu }_3a^3e^{kz}+\mathrm{}]+\mathrm{c}.\mathrm{c}.,$$
(10)
where $`e^{i(kx\omega t)}`$, “c.c.” denotes the complex conjugate, and $`\mu _2`$, $`\mu _3`$, $`\nu _1`$, $`\nu _2`$, $`\nu _3`$ and $`\stackrel{~}{\nu }_3`$ are constants to be determined. We can take $`a`$ to be real without loss of generality. $`\omega `$ must depend on $`a`$ if the appearance of secular terms at the third order is to be avoided. Again expanding in powers of $`a`$, we have

$$\omega =\omega _0+a^2\omega _2+\mathrm{},$$
(11)
where $`\omega _0=\beta E_0k^2`$ is the linear dispersion relation for $`h_0\mathrm{}`$ and $`\omega _2=\omega _2(k)`$.

Our expansion for $`\varphi `$ satisfies the Laplace equation. Since we are seeking a periodic wave train, the boundary conditions (4) and (5) do not apply. Inserting our expansions (9), (10) and (11) into Eqs. (2) and (6), setting $`\alpha =12c`$, and working to third order in $`a`$, we obtain

$$\omega _2=\beta E_0\left(1+\frac{3}{2}\alpha \alpha ^2\right)k^4,$$
(12)
$$\mu _2=\frac{1}{2}\alpha k,$$
(13)
$$\mu _3=\frac{1}{4}\alpha (2\alpha 1)k^2,$$
(14)
$$\nu _1=iE_0,$$
(15)
$$\nu _2=i\left(1\frac{1}{2}\alpha \right)E_0k,$$
(16)
$$\nu _3=\frac{1}{2}i\left(3\frac{7}{2}\alpha +\alpha ^2\right)E_0k^2,$$
(17)
and

$$\stackrel{~}{\nu }_3=\frac{1}{2}i(3\alpha 1)E_0k^2.$$
(18)
The nonlinear dispersion relation is

$$\omega =\beta E_0k^2\left[1\left(1+\frac{3}{2}\alpha \alpha ^2\right)k^2a^2+\mathrm{}\right].$$
(19)
For the isotropic case $`\alpha =1`$, the magnitude of the phase velocity

$$v_p\frac{\omega }{k}=\beta E_0k\left(1\frac{3}{2}k^2a^2+\mathrm{}\right)$$
(20)
is a decreasing function of $`a`$ when $`a`$ is small. In contrast, the phase velocity of small amplitude gravity waves on deep water increases with amplitude. For a given amplitude $`2a`$, the phase velocity is smallest for $`\alpha =3/4`$, or, equivalently, for $`c=1/8`$.

The form of the periodic wave train is given by

$$\zeta =2a\mathrm{cos}(kx\omega t)+\alpha ka^2\mathrm{cos}2(kx\omega t)+\frac{1}{2}\alpha (2\alpha 1)k^2a^3\mathrm{cos}3(kx\omega t)+\mathrm{}$$
(21)
To third order in $`a`$, the wave is symmetric about the crests, even though the applied electric field breaks the right-left symmetry.

IV. Slow Modulation of Finite Amplitude Periodic Wave Trains

We will now relax the simplifying assumptions of Section III. Suppose $`h_0`$ is finite and that the amplitude is a slowly varying function of position. We set

$$\zeta =\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\zeta _nE^n$$
(22)
and

$$\varphi =\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\varphi _nE^n,$$
(23)
where $`Ee^{i(kx\omega _0t)}`$ and $`\omega _0`$ is given by the linear dispersion relation for films of finite thickness, Eq. (7). Since $`\zeta `$ and $`\varphi `$ are real, $`\zeta _n=\zeta _n^{}`$ and $`\varphi _n=\varphi _n^{}`$ for all $`n`$.

Because the amplitude of the wave is of order $`ϵ`$, the coefficients $`\zeta _n`$ and $`\varphi _n`$ are of order $`ϵ^n`$. We may therefore write

$$\zeta _n=\underset{j=n}{\overset{\mathrm{}}{}}ϵ^j\zeta _{nj}$$
(24)
and

$$\varphi _n=\underset{j=n}{\overset{\mathrm{}}{}}ϵ^j\varphi _{nj}.$$
(25)
Here $`n0`$ and the $`\zeta _{nj}`$’s and $`\varphi _{nj}`$’s do not depend on $`ϵ`$. Note as well that $`\zeta _{00}=\varphi _{00}=0`$ and that the $`\zeta _{0j}`$’s and $`\varphi _{0j}`$’s are all real.

To find a steady, periodic wave train (the analog of the Stokes wave), we would take the amplitudes $`\zeta _{nj}`$ to be constants and the $`\varphi _{nj}`$’s to depend on $`z`$ alone. We will go further, and allow the $`\zeta _{nj}`$’s to vary slowly in space and time as viewed from the frame of reference moving with the group velocity. As a result, we will be able to find the steady, periodic wave train and investigate its stability as well. (Our analysis closely parallels Davey and Stewartson’s treatment of modulated gravity waves DS .) To be explicit, we set

$$\zeta _{nj}=\zeta _{nj}(\xi ,\tau )$$
(26)
and

$$\varphi _{nj}=\varphi _{nj}(\xi ,z,\tau ),$$
(27)
where $`\xi ϵ(xv_gt)`$ and $`\tau ϵ^2t`$. The method of multiple scales will be employed, i.e., $`x`$, $`\xi `$, $`t`$, and $`\tau `$ will be treated as independent variables. Using the chain rule, we see that

$$\frac{d}{dx}=\frac{}{x}+ϵ\frac{}{\xi }$$
(28)
and

$$\frac{d}{dt}=\frac{}{t}ϵv_g\frac{}{\xi }+ϵ^2\frac{}{\tau }.$$
(29)
We now insert our expansions for $`d/dx`$ and $`\varphi `$ into the Laplace equation (1). This results in a series of ordinary differential equations for the $`\varphi _{nj}`$’s. Applying the boundary condition (3), we find that

$$\varphi _{11}=A\frac{\mathrm{cosh}kz^{}}{\mathrm{cosh}kh_0},$$
(30)
$$\varphi _{22}=F\frac{\mathrm{cosh}2kz^{}}{\mathrm{cosh}2kh_0},$$
(31)
$$\varphi _{12}=G\frac{\mathrm{cosh}kz^{}}{\mathrm{cosh}kh_0}iA_\xi \frac{z^{}\mathrm{sinh}kz^{}}{\mathrm{cosh}kh_0},$$
(32)
and

$$\varphi _{13}=\frac{1}{2}A_{\xi \xi }(z^{})^2\frac{\mathrm{cosh}kz^{}}{\mathrm{cosh}kh_0}iG_\xi z^{}\frac{\mathrm{sinh}kz^{}}{\mathrm{cosh}kh_0}+H\frac{\mathrm{cosh}kz^{}}{\mathrm{cosh}kh_0},$$
(33)
where $`A`$, $`F`$, $`G`$ and $`H`$ depend only on $`\xi `$ and $`\tau `$. It also follows that $`\varphi _{01}`$ and $`\varphi _{02}`$ are functions of $`\xi `$ and $`\tau `$ alone, while

$$\frac{\varphi _{03}}{z}=z^{}\frac{^2\varphi _{01}}{\xi ^2}.$$
(34)
The next step is quite laborious but is nonetheless straightforward, and so I will omit the details. Inserting our expansions for $`d/dx`$, $`d/dt`$, $`\zeta `$ and $`\varphi `$ into Eqs. (2) and (6) and working to second order in $`ϵ`$, we find that

$$\zeta _{01}=0,$$
(35)
$$\zeta _{11}=i\frac{\sigma }{E_0}A,$$
(36)
$$\zeta _{12}=i\frac{\sigma }{E_0}G+\frac{h_0}{E_0}A_\xi ,$$
(37)
$$\zeta _{22}=\frac{k}{\sigma E_0^2}\left[1\left(1\frac{1}{2}\alpha \right)\sigma ^2\right]A^2$$
(38)
and

$$F=\frac{ik}{2E_0\sigma ^2}(1+\sigma ^2)\left[1\left(2\frac{1}{2}\alpha \right)\sigma ^2\right]A^2.$$
(39)
Equating the coefficient of $`ϵ^3E^0`$ to zero in Eqs. (2) and (6) yields two linear equations for $`\zeta _{02}`$ and $`\varphi _{01}/\xi `$. Solving these for $`\zeta _{02}`$, we obtain

$$\zeta _{02}=\left(\frac{\sigma k}{E_0^2}\right)\frac{\alpha \sigma h_0k2}{v_g/v_01}|A|^2.$$
(40)
When we equate the coefficient of $`ϵ^3E^1`$ to zero in Eqs. (2) and (6), we obtain two equations that are consistent only if

$$iA_\tau =\mathrm{\Gamma }A_{\xi \xi }+\mathrm{\Lambda }|A|^2A,$$
(41)
where the constants $`\mathrm{\Gamma }`$ and $`\mathrm{\Lambda }`$ are given by

$$\mathrm{\Gamma }=\frac{1}{2}\frac{d^2\omega _0}{dk^2}=\frac{\beta E_0}{\sigma }\left(1+2\sigma h_0k\frac{v_g}{v_0}\right)$$
(42)
and

$$\mathrm{\Lambda }=\frac{\beta k^4}{\sigma E_0}\left[\frac{2}{\sigma ^2}1\alpha +\left(\frac{5}{2}\alpha \alpha ^2\right)\sigma ^2\left(1\sigma ^2+\frac{\sigma }{h_0k}\right)\frac{2v_g/v_0\alpha \sigma h_0k}{v_g/v_01}\right].$$
(43)
Eq. (41) is the nonlinear Schrödinger equation, and describes the time development of the amplitude modulation. It has nonlinear plane wave solutions

$$A=A_0e^{i(\mathrm{\Omega }\tau \kappa \xi )},$$
(44)
where $`A_0`$ and $`\kappa `$ are constants and $`\mathrm{\Omega }=\mathrm{\Gamma }\kappa ^2\mathrm{\Lambda }|A_0|^2`$. Let us recast the solution with $`\kappa =0`$ in terms of the original, physical variables. For $`\kappa =0`$, Eq. (37) has the solution $`\zeta _{12}=G=0`$, and we will specialize to this case. We set $`\sigma ϵA_0=iE_0a`$, where $`a`$ is a real constant of order $`ϵ`$. To second order in $`a`$,

$$\zeta =2a\mathrm{cos}(kx\omega t)+\left(\frac{k}{\sigma }\right)\frac{\alpha \sigma h_0k2}{v_g/v_01}a^2+\frac{2k}{\sigma ^3}\left[1\left(1\frac{1}{2}\alpha \right)\sigma ^2\right]a^2\mathrm{cos}2(kx\omega t),$$
(45)
where

$$\omega =\omega _0+\mathrm{\Lambda }\left(\frac{E_0}{\sigma }\right)^2a^2.$$
(46)
Eq. (45) describes a periodic wave train of finite amplitude propagating over a metal or semiconducting layer of finite thickness. Eq. (46) is the nonlinear dispersion relation.

In the limit $`h_0k1`$, Eq. (45) becomes

$$\zeta =2a\mathrm{cos}(kx\omega t)+\alpha ka^2\mathrm{cos}2(kx\omega t)+\frac{1}{2}\alpha ka^2+\mathrm{}$$
(47)
and Eq. (46) reduces to Eq. (19). Eq. (47) reduces to Eq. (21) after a trivial vertical shift of the origin, and so the results of this section have the correct thick film limit.

For $`h_0k1`$, Eq. (45) is

$$\zeta =2a\mathrm{cos}(kx\omega t)\frac{2}{h_0^3k^2}a^2\left[1\mathrm{cos}2(kx\omega t)\right]+\mathrm{}$$
(48)
where

$$\omega =v_0k\left(1+\frac{1}{3}h_0^2k^2+\frac{2}{h_0^4k^2}a^2+\mathrm{}\right).$$
(49)
Eqs. (48) and (49) are valid for $`ka(kh_0)^31`$. It has been shown elsewhere soliton  that the equation of motion for small amplitude, long waves is the Korteweg-de Vries (KdV) equation

$$\zeta _t=v_0\zeta _x+\frac{1}{3}v_0h_0^2\zeta _{xxx}+2\frac{v_0}{h_0}\zeta \zeta _x.$$
(50)
The KdV equation has finite amplitude, periodic solutions (cnoidal waves) of the form

$$\zeta =2a+4a\mathrm{cn}^2[\pi ^1K(m)(kx\omega t)|m],$$
(51)
where cn is a Jacobian elliptic function, $`K(m)`$ is the complete elliptic integral of the first kind,

$$\omega =v_0k\left[1+\frac{4a}{3h_0}\left(\frac{2}{m}1\right)\right],$$
(52)
and $`m`$ is given implicitly by the following relation:

$$mK^2(m)=\frac{2\pi ^2a}{h_0^3k^2}.$$
(53)
Eqs. (51) - (53) give a good approximate solution to the original equations of motion in the limit that $`a/h_0`$ and $`h_0k`$ tend to zero with $`a/h_0^3k^2`$ remaining finite as the limit is taken. If the amplitude of the cnoidal wave is small enough that $`ka(kh_0)^31`$, Eqs. (51) and (52) reduce to Eqs. (48) and (49). Thus, in the thin film limit, the nonlinear wave (45) is a cnoidal wave of small amplitude.

The periodic wave train (44) is marginally stable if $`\mathrm{\Gamma }\mathrm{\Lambda }0`$ but is unstable if $`\mathrm{\Gamma }\mathrm{\Lambda }>0`$ HO . This means that for a given value of $`\alpha `$, the periodic wave train is unstable if $`kh_0`$ exceeds a critical value we shall denote by $`\chi _c=\chi _c(\alpha )`$, and is marginally stable for $`kh_0\chi _c`$. The critical value $`\chi _c=1.609927\mathrm{}`$ for the isotropic case ($`\alpha =1`$) and is infinite for $`\alpha 2`$ and $`\alpha 1/2`$. Fig. 2 is a plot of $`1/\chi _c(\alpha )`$ as a function of $`\alpha `$. It is intereresting to note that the smallest value of $`\chi _c`$ is obtained for a value of $`\alpha `$ differing just slightly from unity ($`\alpha =0.972\mathrm{}`$).

For $`kh_01`$, the product $`\mathrm{\Gamma }\mathrm{\Lambda }`$ is negative, regardless of the value of $`\alpha `$. The periodic wave train is therefore marginally stable in this limit. This is in accord with the fact that cnoidal waves are marginally stable Drazin .

IV. Conclusions

In this paper, the electromigration-induced dynamics of small-amplitude disturbances on a low index crystal surface were studied. Periodic wave trains of finite amplitude were found, as well as their dispersion relation. These wave trains are unstable for $`kh_0>\chi _c(\alpha )`$, and are otherwise marginally stable. The value of the parameter $`\alpha `$ depends on strength of the material anisotropy, and is equal to 1 if the anisotropy is negligible. For $`\alpha =1`$, $`\chi _c=1.609927\mathrm{}`$. For $`\alpha 1/2`$ and $`\alpha 2`$, on the other hand, $`\chi _c`$ is infinite and the wave train is marginally stable no matter what its wavelength.

The equation of motion for slow modulations of the finite amplitude, periodic wave train was shown to be the nonlinear Schrödinger equation. As is well known, this equation is exactly solvable for initial conditions that vanish with sufficient speed as $`|x|\mathrm{}`$ ZS . This exact solution has an number of consequences for electromigration-induced propagation of surface waves. If $`kh_0>\chi _c(\alpha )`$, a localized initial wave packet of arbitrary shape will eventually disintegrate into a number of envelope solitons and a small oscillatory tail. These solitons survive collisions with each other with no permanent change except a shift in phase and position.

What are the prospects of an experimental test of the theory developed in this paper? Quite recently, the SEM-induced dynamics of metal surfaces of large area have received some attention, but unfortunately these experiments were carried out on polycrystalline films, and so they do not permit a test of the theory Shimoni ; Shimoni2 .

The experimental situation is more promising in the case of semiconductors. A number of beautiful experimental studies of the SEM-induced dynamics of single crystal silicon films have been carried out, but on initially planar, vicinal surfaces Latyshev ; bunching . These studies revealed that a vicinal surface is unstable against step-bunching for one current direction, and that it is stable for the opposite current direction. Experimental studies of the dynamics of perturbed, low-index silicon surfaces that are subject to high electric fields have not yet been carried out.

To test the theory, a singly-periodic periodic wave structure could be etched into a low-index silicon surface. After perturbing the surface in this way, a high electrical current would be applied parallel to the ripple wavevector. The subsequent surface dynamics could be imaged by scanning tunneling microscopy, as has already been done for vicinal surfaces.

The analog of the Stokes wave would be obtained by etching many parallel ripples into the wafer surface before the application of current. In this way, the amplitude dependence of the phase velocity could be determined and compared with the theoretical prediction. To observe the development of an envelope soliton (or solitons), on the other hand, a few tens or hundreds of ripples would be etched into the film initially. This would produce a “wave packet” that would ultimately disintegrate into one or more envelope solitons and a small oscillatory tail under the action of an applied current.


text:
# 1 Introduction

## 1 Introduction

Most of the unknown parameters of the Standard Model come from fermion masses and mixing angles. They are put in by hand, and we have no real understanding of their origins. Many models that try to explain their values exist, most involving additional symmetries, but more experimental data will be needed to distinguish between them. Much of the excitement concerning neutrino masses and mixing angles is caused by the hope that their values will facilitate our understanding of the flavor problem.

A potentially valuable source of information about the flavor problem may come from electric and magnetic dipole moments. Just as the Yukawa couplings form a matrix in generation space, the interaction of two fermions with a photon will also be a matrix in generation space. The real and imaginary parts of the diagonal elements will lead to the magnetic and electric dipole moments; the off-diagonal elements wil lead to flavor-changing radiative decays. If one can measure these moments, one would obtain valuable hints of physics beyond the Standard Model and would learn more about the origin of flavor.

This paper is concerned with the electric dipole moments (EDMs) of leptons. There is a possibility of substantial improvement in the experimental bounds on EDMs in the near future. Proposals exist to lower the current bound on the muon’s EDM by six orders of magnitude, and substantial improvement in the bound on the electron EDM is possible. Combining limits on the weak dipole moment of the tau with $`U(1)`$ symmetry improves the current bound by a factor of thirty.

How big might one expect the EDMs to be? In the Standard Model, they are negligibly small. However, they can be much larger in extensions of the Standard Model. In multi-Higgs models, the EDM of the muon can easily be as large as $`10^{24}`$ e-cm, within reach of planned experiments. In leptoquark models, the muon and tau EDM’s are typically $`10^{24}`$ e-cm and $`10^{19}`$ e-cm, respectively. Left-right models have a muon EDM which is typically of the order of $`10^{24}`$ e-cm, and in the minimal supersymmetric standard model(MSSM), the electron EDM is above the experimental bounds if the phases are all unity. The point is that a wide variety of models give EDMs that can be observed in the near future.

Babu, Barr and Dorsner have discussed how the EDMs of leptons scale with the lepton masses. In many models, such as the Standard Model and the MSSM, they scale linearly with the mass. However, in a large number of models, such as multi-Higgs, leptoquark and some flavor symmetry models, the EDM scales as the cube of the lepton mass. In these models, the tau EDM will be 5000 times larger than that of the muon (it should be noted that the electron EDM in some of these models receives a two-loop contribution which varies only linearly, and thus it need not be negligible).

Given cubic scaling, the EDM of a heavy lepton (charged or neutral) could be quite large. For example, a $`100`$ GeV heavy lepton would have an EDM a billion times larger than the muon’s. If the latter is in the expected range of $`10^{24}`$ e-cm, then such a heavy lepton would have an EDM as large as $`0.01`$ e-fermi. This EDM would then dominate the electromagnetic interactions of these leptons, drastically changing their phenomenology. In particular, if a heavy neutrino acquired such an EDM, it could even leave an observable ionization track.

Although originally motivated by cubic scaling models, the possibility of a heavy lepton having a huge EDM is worth studying in a model-independent way. We will see below that plausible models exist in which a heavy neutrino has a very large EDM, even though these models do not have cubic scaling. In this Article, we will concentrate on the possibility that a heavy neutrino could have a large EDM. In an earlier paper by one of us, the differential production cross section for heavy leptons (charged and neutral) with large EDMs was studied, and the fact that heavy neutrinos with large EDMs could leave an ionization track was discussed (although no explicit models were mentioned). This Article extends this earlier work substantially. In Section II, we explicitly present a model in which the EDM of a heavy neutrino can be of $`O(10^{16})`$ e-cm. This will serve as an existence proof that relatively simple models can exist in which this occurs. Current bounds on such large EDMs will be noted. In Section III, the cross section for heavy neutrino production is presented, and the Bethe-Block formula for a large-EDM heavy neutrino traversing matter will be calculated. This extends slightly previous work. The result will be that an ionization track is observable, but one would need to either modify existing detectors or construct a new detector. Since it is unlikely that this would occur without some preliminary evidence that such a large-EDM neutrino exists, we look, in Section IV, at the process of initial and final state radiation, $`e^+e^{}N\overline{N}\gamma `$, which would give an enormous rate for single-photon plus missing energy events, and would provide impetus for looking explicitly for large-EDM heavy neutrinos. Finally, in Section V, we discuss whether these heavy neutrinos could be candidates for the ultra-high energy cosmic rays, and also consider the implications for baryogenesis through leptogenesis, and in Section VI, we present our conclusions.

## 2 Plausible Models

Are there plausible models in which a huge EDM for a heavy neutrino occurs? If one does not rely on any particular model, then a simple effective field theory argument can be used. Suppose one assumes that $`CP`$ violation is due to some sort of new physics at the TeV scale. Then one can write the effective low-energy, dimension-five Lagrangian as

$$=\frac{c}{\mathrm{\Lambda }}\overline{L}_L\sigma _{\mu \nu }i\gamma _5L_RF_{\mu \nu }.$$
(1)
If $`\mathrm{\Lambda }`$ is $`O(1)`$ TeV, and the unknown coefficient is of $`O(1)`$, then this yields a very large EDM of approximately $`10^{15}`$ e-cm for the lepton.

Of course, in realistic models, one expects the EDM to occur via a loop, suppressing the unknown coefficient, and there may not even be a dimension-five operator. What about specific models currently in the literature? A model with cubic scaling is that of Bernreuther, Schroder and Pham, in which CP violation occurs in the Higgs sector. In this model, the EDM of a heavy lepton (they consider the top quark, but the results are unchanged) is constrained by the electron EDM, and the maximum value is $`O(10^{17})`$ e-cm. In other models, such as the model of Babu, Barr and Dorsner in which their parameter $`c=0`$, the electron EDM does not give any such constraints and a large EDM is allowed.

Here, we will provide an explicit example of a model with a large EDM. Consider a model containing singlet fermions, $`E_{L,R}`$ and $`N_{L,R}`$ and two charged singlet scalar fields, $`h_i`$. The fermions can couple to the charged scalars, $`h_i,i`$=1, 2, in the following form

$$Y_i\overline{E}_RN_Lh_i+Z_i\overline{E}_LN_Rh_i+h.c.$$
(2)
where $`Y_i`$ and $`Z_i`$ are complex Yukawa couplings and can generally include flavor indices. The $`h_i`$ can also mix with each other and can be diagonalized by a complex unitary matrix U: $`h_i=U_{ij}H_j`$, with mass eigenstates $`H_i`$. Two such charged singlets are necessary so that a phase redefinition cannot eliminate the CP violation. For simplicity in presentation, we will look at the single $`H`$ case, and just assume that the phase redefinition does not occur. In this case, the coupling is

$$\overline{E}(SP\gamma _5)NH+h.c.$$
(3)
where S and P are arbitrary complex couplings.

It is straightforward to calculate the EDM due to the above interaction. The one-loop feynman diagrams are listed in Fig. 1. The EDM is defined at $`q^2=0`$, and is given by

$$d_N(q^2=0)=e\frac{m_E}{8\pi ^2}Im(SP^{})_0^1𝑑z\frac{(1z)(1+2z)}{z(1z)m_N^2+(1z)m_E^2+zm_H^2}$$
(4)
which, numerically, in units of e-cm, is

$$3.2\times 10^{17}\left(\frac{100\mathrm{G}\mathrm{e}\mathrm{V}}{M_E}\right)\frac{Im(SP^{})}{4\pi }_0^1𝑑z\frac{(1z)(1+2z)}{z(1z)a^2+(1z)+zb^2}$$
(5)
where $`aM_N/M_E`$ and $`bM_H/M_E`$.

The integral is plotted in Fig. 2. We see that EDM’s of $`10^{16}`$ e-cm are quite possible, especially in the region of paramter-space in which the $`E`$ is heavier than the $`N`$ or $`H`$, if one chooses $`\frac{Im(SP^{})}{4\pi }1`$, and can even be an order of magnitude higher if one allows for some fine-tuning. We have not considered the unlikely region of parameter space in which $`M_N>M_E+M_H`$, since the width of the $`N`$ in this case is comparable to its mass, and the definition of the EDM is no longer straightforward.

We thus have presented a model in which the EDM of a heavy neutrino is very large, of $`O(10^{16})`$ e-cm. While this model is, admittedly, designed to have such a large EDM, its ingredients are very simple, and one can imagine that it may be part of a more complicated model. The main point is that having a heavy neutrino with a large EDM is certainly not excluded, and in this paper we are exploring the phenomenological implications.

Are there any current phenomenological bounds? If the heavy neutrino and charged lepton are part of a fourth chiral fermion family (albeit with a right-handed neutrino), then the S and T parameters will cause severe constraints, but if they are part of a vectorlike multiplet, there will be no such constraints. Very explicit models may be constrained, but there are no general bounds. One possible concern would be g-2 of the muon. Here, there will be an effect of a large EDM on the photon propagator. At the lowest order, the feynman diagram is illustrated in Fig. 3. As in the calculation of hadronic contributions to $`(g2)_\mu `$, we calcualte the effect of a large EDM on $`(g2)_\mu `$ by the dispersion integral

$$a_\mu ^{EDM}=\frac{1}{4\pi ^3}_{4m_N^2}𝑑sK(s)\sigma (s)_{e^+e^{}NN},$$
(6)
where K(s) is given in Ref.  and $`\sigma (s)_{e^+e^{}NN}`$ is given in Ref. , as well as below. We cutoff the upper bound of the integral. Letting $`m_N=100`$ GeV, $`a_\mu ^{EDM}=5.9\times 10^{11}(1.3\times 10^{10})`$ for a cutoff of $`1(10)`$ TeV. This is much smaller than the theoretical and experimental uncertainties, and is thus negligible.

## 3 Direct detection

The most dramatic effect of a large EDM of a heavy neutrino will be in the production cross-section and angular distribution. Escribano and Masso noted that the relevant $`U(1)`$ invariant effective operator is given by $`\overline{L}_L\sigma ^{\mu \nu }i\gamma _5L_RB_{\mu \nu }`$, where $`B_{\mu \nu }`$ is the $`U(1)`$ field tensor. This gives a coupling to the photon, which we define to be the EDM, as well as a coupling to the $`Z`$ which is the EDM times $`\mathrm{tan}\theta _W`$. We will include this coupling to the $`Z`$. It turns out that the contribution from the $`Z`$ has very little effect on the numerical results. In general, one should also include an operator coupling to the $`SU(2)`$ field tensor, leading to a very different value for the CP-violating coupling to the $`Z`$. Rather than deal with two parameters, however, we just assume that the latter operator is smaller. If this assumption is false, then it will just make the cross-section even bigger, unless there is fine-tuning. Also, unless the CP-violating coupling to the $`Z`$ is surprisingly large, it will have very little effect on the results.

A discussion of the differential cross section for a heavy charged lepton can be found in Ref. . Here we are interested in heavy neutrino production. The differential cross-section is given by

$$\frac{d\sigma }{d\mathrm{\Omega }}=\frac{\alpha ^2}{4s}\sqrt{1\frac{4M^2}{s}}\left(A_1+\frac{1}{8\mathrm{sin}^42\theta _W}P_{ZZ}A_2+\frac{(14\mathrm{sin}^2\theta _W)\mathrm{tan}\theta _W}{\mathrm{sin}^22\theta _W}P_{\gamma Z}A_3\right)$$
(7)
where

$`A_1`$ $`=`$ $`D^2s\mathrm{sin}^2\theta (1+{\displaystyle \frac{4M^2}{s}})`$
$`A_2`$ $`=`$ $`1+\mathrm{cos}^2\theta {\displaystyle \frac{4M^2}{s}}\mathrm{sin}^2\theta +8C_V\mathrm{cos}\theta +D^2s\mathrm{tan}^2\theta _W(\mathrm{sin}^2\theta +{\displaystyle \frac{4M^2}{s}}(1+\mathrm{cos}^2\theta ))`$
$`A_3`$ $`=`$ $`4D^2s(\mathrm{sin}^2\theta +{\displaystyle \frac{4M^2}{s}}(1+\mathrm{cos}^2\theta ))`$ (8)
where we have dropped the numerically negligible $`C_V^2`$ terms, for simplicity.

The differential and total cross-sections are given in Figures 4 and 5, for a heavy neutrino mass of $`100`$ GeV. We see that for $`D=0`$, the usual $`1+cos^2\theta +C\mathrm{cos}\theta `$ distribution for a lepton (where the $`\mathrm{cos}\theta `$ term is due to $`\gamma Z`$ interference) is found. For a very large EDM, the distribution is completely dominated by the $`\mathrm{sin}^2\theta `$ contribution from the electric dipole term. The angular distribution begins to deviate significantly from the $`D=0`$ distribution for $`D`$ greater than $`3\times 10^{17}`$, and by $`D=3\times 10^{16}`$ the distribution is very close to $`\mathrm{sin}^2\theta `$. One can see that EDM’s of the size noted earlier will dramatically alter the angular distribution.

Note that since the cross-section varies as $`D^2`$, for an EDM as ridiculously large as $`1.0`$ e-fermi, it would be almost a microbarn!! Cross sections this large will violate unitarity. The unitarity limit can be estimated by setting the relevant effective interaction strength, $`\alpha D\sqrt{s}`$, equal to unity. For $`\sqrt{s}=500`$ GeV, $`D=10^{15}`$ e-cm, and $`\alpha 1/125`$, this effective interaction strength is $`0.25`$. Thus the unitarity limit will not be reached unless the EDM is larger than $`10^{15}`$ e-cm, and so the behavior discussed above will appear for a wide range of parameter space without violating the unitarity bound. Another way of saying this is to note that the larger the EDM, the smaller the scale at which the physics responsible for the effective interaction sets in, and for an EDM larger than $`10^{15}`$ e-cm, that scale is less than $`\sqrt{s}`$.

Production of a fermion with a large EDM has been considered elsewhere. It is given in the context of a top quark EDM by Bernreuther et al., and given in the context of tau-pair production in Ref. . This latter paper noted how one can use CP-odd angular correlations to search for a tau EDM, and this method has been used by experimentalists. However, there has not been any discussion of the possibility of a large EDM for heavy neutrinos, and here we see a unique and unusual signature.

What is the purpose of this calculation? After all, it would not seem remotely possible to detect heavy neutrinos directly. This calculation would only be meaningful if the heavy neutrino were heavier than the charged lepton, and could thus be detected via its decay. Here, one could again look for CP-odd correlations, as discussed in Ref. . However, most models have the charged lepton heavier than the neutrino, and thus the decay can only occur through mixing with the very light neutrinos. As discussed in detail in Ref. , this mixing could be very small, and these heavy neutrinos could be effectively stable. As we will see below, however, it may be possible to detect these neutrinos initially through final state radiation, and then directly.

Because the heavy neutrino has such a large electric dipole moment, it might be possible to detect it directly through its ionization loss. To estimate whether or not this might be possible, we follow the derivation of the classical Bohr formula outlined in Jackson, replacing the electric field from a charge with the electric field from a dipole. This is a classical calculation; if one is not too far above threshold, the neutrino will be moving non-relativistically. Consider a heavy neutrino moving in the $`x`$-direction, and an electron at an impact parameter $`y=b`$. The impulse given to the electron is $`\mathrm{\Delta }\stackrel{}{p}=_{\mathrm{}}^{\mathrm{}}e\stackrel{}{E}𝑑t`$. The result will depend on the orientation of the dipole. Suppose that the dipole is in the $`z`$-direction, transverse to the plane of the particle motion and the electron. Then the only nonzero electric field component is in the $`z`$-direction, and $`E_z=\frac{eD}{4\pi ϵ_o}(b^2+v^2t^2)^{3/2}`$, where the time of closest approach is defined to be $`t=0`$, $`v`$ is the velocity, and $`eD`$ is the size of the electric dipole moment. Integrating, one finds that the impulse is $`\frac{e^2D}{4\pi ϵ_o}\frac{2}{vb^2}`$. Suppose that the dipole is in the $`y`$-direction. Then the electric field components are $`E_x=\frac{eD}{4\pi ϵ_or^3}(3\mathrm{sin}\theta \mathrm{cos}\theta )`$ and $`E_y=\frac{eD}{4\pi ϵ_or^3}(3\mathrm{cos}^2\theta 2)`$, where $`r^2=b^2+v^2t^2`$ and $`\mathrm{tan}\theta =b/vt`$. Integrating, the impulse in the $`x`$-direction vanishes, as expected by symmetry, and the impulse in the $`y`$-direction is also $`\frac{e^2D}{4\pi ϵ_o}\frac{2}{vb^2}`$. Finally, if the dipole is in the $`x`$-direction, the electric field components are $`E_y=\frac{eD}{4\pi ϵ_or^3}(3\mathrm{sin}\theta \mathrm{cos}\theta )`$ and $`E_x=\frac{eD}{4\pi ϵ_or^3}(3\mathrm{cos}^2\theta 1)`$. Both of these integrate to zero, so there is no net momentum transfer in this case. Since the orientation is generally arbitrary, and we are only interested in a rough order-of-magnitude estimate, we will take the impulse to be $`|\mathrm{\Delta }\stackrel{}{p}|=\frac{e^2D}{4\pi ϵ_o}\frac{2}{vb^2}`$. Note that the impulse from an electric charge is just $`\frac{2e^2}{4\pi ϵ_ovb}`$, so this result is expected dimensionally. This impulse is then converted into an energy transfer, $`\mathrm{\Delta }E=|\mathrm{\Delta }\stackrel{}{p}|^2/2m`$. Jackson notes that the maximum energy transfer is $`\mathrm{\Delta }E_{max}=2m\gamma ^2v^2`$, and thus the minimum impact parameter is $`b_{min}^2=e^2D/(m\gamma v^2)`$. The energy loss is obtained by cylindrically integrating over the impact parameter

$$\frac{dE}{dx}=2\pi NZ_{b_{min}}^{\mathrm{}}\mathrm{\Delta }E(b)b𝑑b$$
(9)
and we obtain

$$\frac{dE}{dx}=2\pi N_A\left(\frac{e^2}{4\pi ϵ_o}\right)D\gamma \frac{Z}{A}$$
(10)
which is the corresponding formula to the Bohr result. Since we are doing a classical calculation, one can set $`\gamma =1`$. Note that the logarithm in the usual Bethe-Block formula is absent and the electron mass and particle velocity drop out. This is due to the extra power of $`b`$ in the expression for the impulse, which arises on dimensional grounds. Plugging in numbers, one finds that $`\frac{dE}{dx}=10^{12}D\mathrm{MeVg}^1\mathrm{cm}^2`$.

For EDM’s of $`10^{16}`$ e-cm, this gives an energy loss which is approximately $`10^4`$ of the usual energy loss for a charged particle. This is challenging for experimenters, but not impossible, requiring specialized detectors. For example, if the heavy neutrino is produced in the decay of a charged lepton, then it would be produced in coincidence with a real or virtual W, which would eliminate backgrounds. Note that the event rate, for $`D=10^{16}`$ cm, is huge, corresponding to $`100,000`$ events for an integrated luminosity of $`100\mathrm{fb}^1`$. This is several events per hour, and thus a specialized detector located on the outside of the main detectors could suffice, without needing a full $`4\pi `$ coverage (note that the $`N`$’s will emerge back-to-back). More importantly, however, is that initial and final state radiation will be huge, and thus standard searches for single-photon plus missing energy events will discover a huge signal, which could then be followed by a specialized search for direct detection. We now turn to the possibility of detection through single photon events.

## 4 Detection through single photon events

In the standard model, the process $`e^+e^{}\nu \overline{\nu }\gamma `$ has been calculated and used to determine the number of light neutrinos . The signature of the process is a single photon from initial state radiation plus missing energy. The process didn’t rule out the existence of heavy neutrinos. Similarly, the process $`e^+e^{}N\overline{N}\gamma `$ can be used to study the effects of large EDM of heavy neutrinos. Now the single photon can also come from final state radiation because of the large EDM.

The interaction of heavy neutrinos with the photon is $`ie\overline{N}D\sigma _{\mu \nu }\gamma _5q^\nu NA^\mu `$, where q is the momentum transfer. The EDM is defined to be $`eD(q^2=0)`$. In this case, the photon momenta at the relevant vertices are different, and are certainly not zero. However, we presume that $`D(q^2)`$ does not vary too rapidly with $`q^2`$, as is the case in which the momentum transfer is not substantially greater than the masses of the particles in the loop. The feynman diagrams are listed in Fig. 6. The contribution due to the exchange of the Z boson is neglected here since it is numerically small.

The calculation is straightforward. It is clear that the rate is the sum of terms proportional to $`D^2`$, $`D^3`$ and $`D^4`$. The contributions proportional to $`D^2`$ , $`D^3`$ and $`D^4`$ are presented separately in Fig. 7 for $`D=10^{16}`$ e-cm. For different $`D`$’s, these curves scale appropriately. One can see that for $`D<6\times 10^{16}`$ e-cm, expected in realistic models, the $`D^2`$ term dominates. Because of this scaling, instead of falling as $`1/s`$, the cross section reaches a constant (it would grow as $`s`$ if the $`D^4`$ term dominated. As discussed above, however, it does not reach the level at which unitarity is a big concern. The total cross section is then given, to a good approximation, by the $`D^2`$ term in Fig. 7.

The cross section will be singular as the $`|\mathrm{cos}\theta _\gamma |`$ tends to 1. We have used cuts in the photon energy and in the angle: $`\frac{E_\gamma }{\sqrt{s}}0.1`$ and $`|\mathrm{cos}\theta _\gamma |0.94`$. For example, at $`\sqrt{s}=500`$ GeV, the cross section is about $`1.5\times 10^4`$ nb (and scales as $`D^2`$). This is a very large event rate, corresponding to about an event every few minutes (for $`D=10^{16}`$ e-cm) at a linear collider with a luminosity of $`3\times 10^{34}\mathrm{cm}^2\mathrm{sec}^1`$.

The angular distribution is also important. It is shown in Fig. 8 for $`D=10^{16}`$ e-cm, and the resulting forward-backward asymmetry is shown as a function of $`D`$ in Fig. 9. The angular distribution of the Standard Model is nearly forward-backward symmetric, but in this case there is a sizable asymmetry. For much of the parameter space, the asymmetry is $`0.48`$.

Clearly, the first signature of an EDM for a heavy neutrino will be a very large event rate for single photon events, with the forward-backward asymmetry measurement giving a confirmation of this particular interpretation of such a large rate.

## 5 Astrophysical and Cosmological Significance

One of the greatest mysteries in cosmic ray astronomy concerns the existence of cosmic rays with energies above the GZK bound. A substantial number of cosmic rays with energies well above $`10^{20}`$ eV have been observed, interacting in the upper atmosphere. Yet a charged particle’s mean free path through cosmic microwave background, at these energies, is only about $`50`$ megaparsecs, and no sources of such high energy particles are known to exist within this distance. Neutral particles will not have this problem. However, the fact that the cosmic rays shower high up in the atmosphere implies that the cross-section must be of strong interaction strength. The only strongly interacting long-lived neutral particle is the neutron, yet its lifetime is too short. Neutrinos will not interact high in the atmosphere, if they have conventional weak interactions.

Several papers have considered the possibility that neutrino interactions might become strong at very high energies, due to exchange of Kaluza-Klein gravitons. One might wonder whether heavy neutrinos with a sufficiently large EDM might interact strongly enough in the upper atmosphere, but not interact strongly with the microwave background. Are these candidates for the ultra-high energy cosmic rays? One can redo the calcuation described above concerning ionization loss in a detector, but now looking at the ultrarelativistic limit. This is a tedious, but straightforward calculation.

In this case, the cross section can not be high enough. Recall that the EDM is defined at $`q^2=0`$. However, here one would expect $`q^2`$ to be enormous. Certainly, in the first few interactions, it will be much larger than the mass-squared of the particles in the loop that produces a large EDM in the first place. As a result, the effective EDM will be much smaller, and the interaction will not be strong. These heavy neutrinos will not interact high in the atmosphere, and may not interact much more than conventional high energy neutrinos.

Are there any other potential astrophysical or cosmological effects? One possibility concerns the generation of the baryon asymmetry. Suppose the heavy neutrino decays radiatively into light neutrinos. Then the huge EDM, which violates CP, may cause the decay rate into neutrinos to be different than the decay rate into antineutrinos. This would result in a large lepton asymmetry. Since non-perturbative sphaleron interactions can convert a lepton asymmetry into a baryon asymmetry, this would result in a baryon asymmetry. The reader is referred to the article of Riotto and Trodden for a nice review. This possibility is currently under investigation.

## 6 Conclusions

In this Article, we have explored the possibility that a heavy neutrino could have a huge electric dipole moment. It has been shown that plausible models exist that would lead to a dipole moment as large as $`10^{16}`$ e-cm, and that this would dramatically alter the electromagnetic properties of these neutrinos. The first signature of such a neutrino would come from a large enhancement of the single photon plus missing energy event rate. Then, one could actually look for the ionization track of these neutrinos in specially designed detectors. Cross sections and angular distributions are calculated as a function of the electric dipole moment. While these neutrinos can not explain the ultra high energy cosmic rays, they may have interesting implications for baryogenesis.

Acknowledgments

We thank Mariana Frank as one of us (Nie) is now doing postdoctoral research under her supervision. We also are extremely grateful to John Ralston, for repeatedly warning us that the EDM at high energies would be softened, making it unlikely that these neutrinos could be the UHE cosmic rays. We finally believed him. We also thank Jack Kossler and Chris Carone for many useful discussions. This work is supported by the National Science Foundation grant number PHY-9900657.