format: text

text:
# On network bipartivity

## I Introduction

Any system, natural or man-made, consisting of entities that interact pairwise can be described in terms of a network. Networks in the real life often contain some degree of randomness, and has also some structure arising from the strategies or laws the entities follow to make new contacts. Such networks—that can only be described as having both randomness and structure—are called complex networks and has lately received much attention in the physicist community review ; review2 . Among the most important developments in this recent surge of activity in network research is arguably the categorization and quantification of static network structures such as clustering WS , degree distribution sf , assortative mixing coefficient assmix , grid coefficient grid , etc. A network with no circuit of odd length is called bipartite. Many systems are naturally modeled as bipartite networks: Biochemical networks can be described by vertices representing chemical substances separated by vertices representing chemical reactions jeong . As another example, we have the so called “two-mode” representation of affiliation networks where one kind of vertices represents e.g. organizations and the other type represents individual actors, and the edges indicates to which organizations an actor belongs. But there are also networks that are not necessarily bipartite, but closer to bipartite than what can be expected from a completely random network. Examples of such networks are those that are formed by two types of agents with a preference for heterophilous interaction (human sexual contacts liljeros ; lea , and human romance or partnership networks partner  being two cases). In many cases one knows the type of the individual vertices (the gender of the actors in the examples above) freeman , but in other cases such information might be lacking (the data studied in Ref. HEL  for a concrete example). Nevertheless, the ‘bipartivity’—how far away from being bipartite a graph is—is a measurable structure; and therefore, we believe, deserves attention.

How can we measure bipartivity? The idea we use in this paper is the following: We suppose that all agents of one type tried their best in forming a connection to an agent of the other type. Then we measure to what extent this assumption fail. We can assign a label $`\sigma _v\{1,+1\}`$ to each vertex $`v`$ and check for the maximal fraction of edges between vertices of different sign. This fraction will be equal to or higher than the actual fraction of edges between vertices of different type. But, at least for strong heterophilous preference in the network formation, the difference should be small. For weak heterophilous preference this approach will likely fail to produce a correct classification of the individual vertices. Still, the number of even circuits should be larger than in a network created under the same circumstances but with no heterophilous preference; and this will (as we will see) give a lower value of such a bipartivity measure. So even if we cannot reproduce the correct fraction of vertices of different type, we have a measure that is a monotonous function of the strength of the heterophilous preference. It is convenient (at least for people familiar with statistical mechanics) to phrase a problem like this in terms of the antiferromagnetic Ising model. Our bipartivity measure—the maximal fraction of edges between vertices of different sign—is directly related to the ground state energy of the antiferromagnetic Ising model (the relation is given in Sect. II.1.1). Throughout the paper we will often use the terminology of such spin systems, such as the antiferromagnetic Ising model. For example we talk of an edge between two vertices of the same tag as a ‘frustrated’ edge.

The spin system analogy to combinatorial optimization problems such as the one we are facing—to find minimal fraction of frustrated edges—is nothing new. With this approach the fraction of frustrated edges defines a cost function corresponding to the energy of the spin system. The two most studied problems in this area are the $`p`$-coloring problem and the graph bisection problem. In the $`p`$-coloring problem the question is whether or not the vertices of a graph can be assigned one of $`p`$ colors in such a way that no edge goes between two vertices of the same color. This problem is solvable in linear time for $`p=2`$, but NP-complete (i.e. in the general case not calculable in polynomial time hope ) for $`p>2`$. The graph bisection problem (also NP-complete) is to partition the vertex-set into two sets of equal size such that the number of edges between the two sets is minimized jerrum ; schreiber ; FA . Both these problems can, just as ours, be phrased in terms of spin-models with antiferromagnetic interaction. Our minimization problem is a little bit different from the bisection problem in that the two sections can have arbitrary sizes. However, as in the bisection and $`p`$-coloring problems we are also faced with an NP-complete optimization problem. (Our aim—to find the ground state energy of antiferromagnetic Ising model can be mapped to a min-flow max-cut problem alava  which is NP-hard on general networks karp .)

As the spin models of statistical physics are familiar to statistical physicists, it is not surprising that topics like the Ising and XY models on various model networks bw ; spin  have received much attention in physicists’ network literature. The motivation for such studies, as models of real-world systems, is that they can capture some features of opinion formation or similar social processes socstatmech . The present work can also be described as a study of a spin model on a complex network, but unlike the above mentioned studies, the spin model is used as a tool to measure a static network structure.

## II The measures

In the following sections we will go through the two bipartivity measures. We state the definitions, dissect the algorithms and give analytic discussions about the limit properties.

We represent a undirected network by $`G=(V,E)`$ and a directed network by $`G_{\mathrm{dir}}=(V,A)`$, where $`V`$ is the set of vertices, $`E`$ is a set of edges (or undirected pairs of vertices), and $`A`$ is a set of arcs (or ordered pairs of vertices). A path of length $`l`$ is a sequence of vertices $`v_1,\mathrm{},v_l`$ such that $`(v_i,v_{i+1})E`$ (or $`(v_i,v_{i+1})A`$ for directed graphs); a circuit is a path where the first and last vertex are identical. In an elementary path, or circuit, no vertex appears twice (except the first and last in case of circuits). In the present paper we will only talk about elementary paths and circuits—so, for brevity we omit the word ‘elementary.’ Throughout the paper, when necessary, we let sub- or superscript ‘dir’ denote directed versions of quantities. In many cases the generalization from undirected to directed networks is straightforward; in these cases we will pursue the discussion in the framework of undirected networks.

### II.1 The measure $`b_1`$

#### II.1.1 Definition

The first measure we consider is simply the fraction of unfrustrated edges in the ground state of the antiferromagnetic Ising model on the network. In terms of the antiferromagnetic Ising model the quantity can be written as

$$b_1=1\frac{M_{\mathrm{fr}}}{M}=\frac{1}{2}\frac{E_0}{2M},$$
(1)
where $`M_{\mathrm{fr}}`$ is the number of frustrated edges in the ground state (the usual cost function in the two-coloring problem). $`E_0`$ is the ground state energy

$$E_0=\underset{\{\sigma _v\}}{\mathrm{min}},H$$
(2)
where $`H`$ is the Hamiltonian of the antiferromagnetic Ising model:

$`H`$ $`=`$ $`{\displaystyle \underset{(v,w)E}{}}\sigma _v\sigma _w`$ (3a)
$`H_{\mathrm{dir}}`$ $`=`$ $`{\displaystyle \underset{(v,w)A}{}}\sigma _v\sigma _w`$ (3b)
The directed quantity is obtained by substituting $`H`$ by $`H_{\mathrm{dir}}`$ in Eqs. (1) and (2), and edges by arcs in the above discussion. The topology of the energy landscape is determined by the underlying network, and can in general be very complex barahona .

#### II.1.2 Limit properties

The $`b_1`$ measure takes values in the interval $`(1/2,1]`$. The upper bound is attained for bipartite graphs. It is easy to see that $`b_1`$ cannot be lower than $`1/2`$: Consider a ground state configuration for which the opposite is true. Then there must be at least one vertex with more than half of its edges frustrated. Flipping this spin would reduce the energy, which contradicts the fact that the system is in the ground state. We do not know if this bound is realized for any finite graphs, but $`b_1=1/2`$ is the limit value for $`b_1`$ for a fully connected graph as $`N\mathrm{}`$: Partition the fully connected graph $`K_N`$ of $`N`$ vertices (and $`M=N(N1)/2`$ edges) into one set of $`N^{}`$ and one set of $`NN^{}`$ vertices and assign opposite spins to the elements of these sets. The number of frustrated edges is precisely the number of edges within each set which is:

$`M_{\mathrm{fr}}(K_N)`$ $`=`$ $`{\displaystyle \frac{N^{}(N^{}1)}{2}}+{\displaystyle \frac{(NN^{})(N1N^{})}{2}}`$ (4)
$`=`$ $`MN^{}(NN^{}).`$
Thus the minimum number of frustrated edges is exactly $`N^2/4N/2`$ for $`N^{}=N/2`$, and the fraction of unfrustrated edges is

$$b_1=\frac{1}{22/N}\frac{1}{2}\text{ as }N\mathrm{}.$$
(5)
The above arguments can be generalized to directed networks straightforwardly.

#### II.1.3 Minimization by exchange Monte Carlo

The complexity of the “energy landscape” of the antiferromagnetic Ising model on an arbitrary network is difficult to judge a priori. There are indications that no natural network would be too hard for a regular simulated annealing approach simann ; jerrum . To be on safer ground, we use a Monte Carlo scheme that is evidently very efficient to sweep even an extremely ‘rugged’ energy landscape without getting stuck in local minima—the so called exchange Monte Carlo (XMC) xmc . The idea of exchange Monte Carlo is to run standard Metropolis Monte Carlo for $`N_T`$ replicas of the system, each at a specific temperature. Then from time to time two replicas at adjacent temperatures are compared, and with a probability

$$P_{\mathrm{exch}.}=\{\begin{array}{cc}1\hfill & \text{if }\mathrm{\Delta }<0\hfill \\ e^\mathrm{\Delta }\hfill & \text{otherwise}\hfill \end{array},$$
(6)
where

$$\mathrm{\Delta }=\left(\frac{1}{T}\frac{1}{T^{}}\right)(E^{}E),$$
(7)
and $`E`$ is the energy of the configuration at temperature $`T`$ (similarly for $`T^{}`$ and $`E^{}`$), and $`T<T^{}`$. the two replicas are swapped between the temperatures. This condition is designed so that the Monte Carlo scheme preserves the Boltzmann distribution. This is not decisive for us who are looking for the ground state energy, rather that performing a proper sampling of the configuration space, but anyway kept in our measurements. Besides just running the XMC scheme we also periodically quench the system, i.e. we sweep through all vertices of the network consecutively and flip spins that lower the energy. The sweeps are continued until a sweep with no spin-flips has occurred. For later reference we introduce the notations $`t_{\mathrm{avg}}`$ for the total number of MC sweeps—we refer to the number of MC sweeps as ‘time’—$`t_{\mathrm{quench}}`$ for the time between each quench, $`t_{\mathrm{exch}}`$ for the time between exchange trials, $`t_{\mathrm{measure}}`$ for the time between measurement sweeps (where the energy is sampled).

For the exchange Monte Carlo scheme to efficiently sample the configuration space all replicas needs to tour the whole range of temperatures in a reasonably short time. At the same time one would not like the exchange trials, at any neighboring temperatures, to be constantly affirmative—then the separation of the two temperatures would be of no use. We follow Ref. xmc  and choose the temperature set

$$T_i=T_{\mathrm{low}}\left(\frac{T_{\mathrm{high}}}{T_{\mathrm{low}}}\right)^{(i1)/(N_T1)},$$
(8)
where $`1iN_T`$ enumerates the replicas. $`T_{\mathrm{low}}`$ is the lowest and $`T_{\mathrm{high}}`$ represent the highest temperatures respectively. To find the actual parameter values (which will be stated in Secs. IV.1 and IV.2) one has to check that the replicas travels throughout the temperature range with reasonable exchange ratios for all temperature gaps.

### II.2 The measure $`b_2`$

Apart from finding an approximative value of $`b_1`$, one can also define a quantity that is exactly solvable in polynomial time. Our intention is in the first hand not to make a heuristic algorithm for calculating $`b_1`$, but rather a quantity that captures the same structure, i.e. that grows monotonously with $`b_1`$.

That a graph contains no odd circuits is the defining property of bipartiteness intro . It is thus natural that we base a bipartivity measure on an odd-circuit count in some way. Unfortunately, defining a quantity in this way becomes a little bit more complicated than at first expected. One complication is that a graph can be very close to bipartite and still contain many odd-circuits (see Fig. 1(a)). A way of dealing with this problem is to mark as few edges as possible such that each odd circuit contains at least one marked edge. In many cases a marked edge will correspond to a frustrated edge of the ground state of the antiferromagnetic Ising model. In Fig. 1(a) only the upper, horizontal edge needs to be marked. Another problem one faces is how to deal with odd circuits of different length—in a network with very few odd circuits a circuit of, say, length seven would contribute as much to the global frustration of the network as a triangle (a subgraph of three adjacent vertices—see Fig. 1(b)). But in many real networks the total length of the odd circuits is very long (this is true for all networks we measure, see Sect. III.2), much larger than $`M`$, in these cases the short circuits are in general the most important in determining the ground state configuration. For example, in Fig. 1(c) $`M=23`$, and while we have 11 triangles, summing the lengths of all odd circuits gives 218 (33 from the 11 triangles, 45 from the nine circuits of length five, and so on). However, only the triangles contributes to the ground state configuration in the sense that each triangle has the same configuration as the ground state of an isolated triangle, while all odd circuits of length larger than four (e.g. the periphery) has not the best coloring for a circulant of that length. To deal with this we need to weight short circuits higher than long. We will do this by assigning a cut-off length and neglect all circuits exceeding this length.

#### II.2.1 Definition

Now, we make an algorithm of the above ideas as follows: Let $`C_n`$ be the set of odd circuits of length $`n`$. Let $`\mathrm{\Sigma }(C_n)`$ be the accumulated length of the circuits in $`C_n`$ (so, for example $`\mathrm{\Sigma }(C_3)=3`$ in Fig. 1(b)). Now we assign the cut-off $`3M`$ to $`\mathrm{\Sigma }(C_n)`$, and let $`\widehat{n}`$ be the smallest $`n`$ such that $`\mathrm{\Sigma }(C_n)3M`$. Next we turn to the marking procedure sketched above. Let $`\nu (e)`$ denote the number of circuits in $`C_{\widehat{n}}`$ passing through the edge $`e`$. Clearly edges of high $`\nu `$ are likely to be frustrated in the ground state (viz. Fig. 1(a)). We now estimate $`M_{\mathrm{fr}}`$ roughly as the number of edges that has to be marked so that each odd circuit of length $`\widehat{n}`$ is marked at least once. To be precise we perform the following algorithm:

1. Start with $`C=C_{\widehat{n}}`$.
2. Sort the edges in order of $`\nu `$.
3. Repeat the following while $`C\mathrm{}`$:

1. Mark the edge $`e`$ with highest $`\nu `$.
    2. Remove all circuits in $`C`$ containing $`e`$.
    3. Recalculate $`\nu `$ for each edge.

Then the number of iterations $`m^{}`$ is the assessment of $`M_{\mathrm{fr}}`$, and we define our bipartivity measure as

$$b_2=1\frac{m^{}}{M}.$$
(9)
This algorithm is not an attempt to actually identify the frustrated edges, rather it is supposed to give a high $`M_{\mathrm{fr}}`$ for a system with high (total) geometric frustration, and vice versa: Firstly, it does not necessarily find the minimal number of edges needed to be marked for all odd circuits of length less than $`\widehat{n}`$ to contain a marked edge. But we expect this steepest descent optimization to come close in most cases. Secondly, an odd circuit can in reality only have an odd number of frustrated edges, but in the algorithm there is no such restriction on the number of marked edges.

In case there are more than one edge with the highest $`\nu `$ (in step 3a of the algorithm) we choose the edge to mark at random. The variance between different random seeds turns out to be negligible in most cases. We will run the algorithm for different seeds to choose the highest $`b_2`$ value, and get an idea about the error in $`b_2`$ from the selection of edge to mark. An alternative (and more ambitious) approach would be to iterate the whole calculation until the highest $`b_2`$ has reappeared a fixed number of times (cf. ww ).

If we assume a sparse network (i.e. $`NM`$) the running time of the algorithm above is $`O(M^2)`$. To see this we first note that there can be at most $`O(M)`$ iterations at step 3. To find the edge with highest $`\nu `$ (in step 3a) we do not need to sort all edges more than once (as done in step 2). Instead we can find this out while recalculating $`\nu `$ (in step 3c). Removing all circuits containing $`e`$ (as in step 3b) can be done in time bounded by the total length of circuits containing $`e`$, which cannot be larger than $`3M`$. Step 3c also needs to go through all circuits passing $`e`$ and thus needs the same running time as step 3b. To sum this up the running time for this section of the algorithm is of order $`N^2`$.

#### II.2.2 Limit properties

In the $`N\mathrm{}`$ limit the $`b_2`$ measure lies in almost the same interval as $`b_1`$. The upper limit $`b_2=1`$ is attained if and only if the graph is bipartite. (If the graph is bipartite $`C_{\widehat{n}}`$ is empty and $`\nu (a)=0`$ for all $`a`$, so $`m^{}=0`$ and $`b_2=1`$. If there exists odd circuits $`m^{}0`$, so $`b_2<1`$.) $`b_2`$ cannot be as low as 0 (if one marks all edges, every circuit must be marked). Since the $`b_2`$-definition is inspired by the ground-state configuration of the antiferromagnetic Ising model, we expect a similar lower bound to $`b_2`$ as to $`b_1`$. In Appendix A we argue that the lower bound on the $`b_2`$, as for the $`b_1`$ measure, is $`1/2`$ in the $`N\mathrm{}`$ limit.

#### II.2.3 The complete algorithm

So far we have overlooked the central part in calculating the $`b_2`$ measure—namely to find odd circuits. To do this we use a modified version of Johnson’s algorithm johnson . In principle Johnson’s algorithm is a depth first search where, to avoid futile searching, some vertices are blocked while stepping down the search tree. The running time for Johnson’s algorithm is $`O(M(C+1))`$ (if $`M>N`$) where $`C`$ is the total number of circuits. Now $`C`$ can grow fast with $`N`$ which would make the finding of all odd circuits a quite intractable computation. In many cases the cut-off of the circuit length, that we introduced above to give less priority to long circuits, saves us by setting a limit on the search depth. To implement this we let $`\overline{n}`$ be the current upper bound on circuit length (or search depth), and $`\overline{\mathrm{\Sigma }}`$ be the current sum of odd circuits $`\overline{n}`$. As soon as $`\overline{\mathrm{\Sigma }}M`$ we iteratively decrease $`\overline{n}`$ by $`2`$ and recalculate $`\overline{\mathrm{\Sigma }}`$ until $`\overline{\mathrm{\Sigma }}<3`$. If $`\overline{\mathrm{\Sigma }}<M`$ when the search is over we rerun the procedure where we use $`\overline{n}+2`$ as our new (fixed) $`\overline{n}`$ note:alt . When the search is over we assign $`\widehat{n}`$ the value $`\overline{n}`$. For dense bipartite graphs the algorithm is intractable. In the worst case, the full bipartite graph, $`K_{N/2,N/2}`$, there are

$$C(K_{N/2,N/2})=\underset{k=4}{\overset{N}{}}\frac{1}{2k}\left[\frac{(N/2)!}{(N/2k/2)!}\right]^2$$
(10)
circuits (where the sum is over even values of $`k`$note:bip  giving a running time of $`O(N^2C(K_{N/2,N/2}))`$. One can of course decide whether or not a graph is bipartite in linear time, but non-bipartite cases of similar complexity are easily constructed (by, e.g., adding an isolated triangle). In practice these worst cases are, probably, very rare—a, relatively speaking, very low density of odd circuits is needed to get a small $`\widehat{n}`$—even in the real-world network with highest bipartivity we have $`\widehat{n}=3`$. In this case ($`\widehat{n}=3`$) all odd circuits are found in $`O(M^2)`$ time.

Now we turn to a more complete description of the algorithm. Johnson’s algorithm takes the ‘least’ (smallest in some enumeration) vertex in a strongly connected subgraph as its starting point. To find strongly connected components we use the algorithm in Ref. SCC . To sum up, the algorithm reads:

1. Mark all vertices as unchecked.
2. While there are unchecked vertices, iterate the following:

1. Pick an unchecked vertex $`v`$.
    2. Find the largest strongly connected component $`\mathrm{\Lambda }_v`$ containing $`v`$.
    3. Set $`\mathrm{\Lambda }:=\mathrm{\Lambda }_v`$ and repeat the following steps as long as $`\mathrm{\Lambda }\mathrm{}`$:

    1. Pick the least vertex $`u`$ of $`\mathrm{\Lambda }`$.
        2. Call a subroutine implementing the modified Johnson’s algorithm. Recalculate $`\overline{n}`$ and add $`C_{\overline{n}}`$ to a list $`𝒞`$. Delete circuits longer than $`\overline{n}`$ from $`𝒞`$.
        3. Delete $`u`$ from $`\mathrm{\Lambda }`$.
3. Set $`\widehat{n}:=\overline{n}`$.
4. Run the algorithm described above (in Sect. II.2.1) to mark edges and calculate $`b_2`$.

In all cases, step 2 sets the limit on running time. As mentioned, in most application we expect the running time of step 2 to be $`O(M^2)`$ (similarly to that of step 4).

## III The Networks

### III.1 Test networks with tunable bipartivity

To test and compare the $`b_1`$ and $`b_2`$ quantities we construct three types of test networks where the bipartivity can be tuned by model parameters. The principle behind all models is to start from bipartite networks and add lesser or greater number of edges within a partition to create odd circuits.

One type (Model 1) is a quite straightforward generalization of the Erdös-Renyi (ER) model ER : We partition the vertices in two disjoint sets of sizes $`\stackrel{~}{N}`$ and $`N\stackrel{~}{N}`$. Then we add $`r_1M`$ edges randomly between vertices of the different sets, and $`(1r_1)M`$ edges regardless of what set the vertices belongs to (see Fig. 2(a)). In this way we interpret $`r_1`$ as the strength of the heterophilous preference in a model where bipartivity is the only structural bias. The choice of vertex pairs is done with randomness, the only restriction being that loops and multiple edges are not allowed. If $`r_1=0`$ the model reduces to the ER model, while for $`r_1=1`$ the networks are bipartite (cf. Ref. nws ). This model is probably the most random (i.e. having least structural biases) model with tunable bipartivity. The disadvantage is that the expectation values of $`b_1`$ and $`b_2`$ are hard to calculate (even in the frustrated limit $`r_1=0`$).

Model 2 interpolates between two-dimensional square- and triangular lattices. We start, for $`r_2=0`$, with a triangular grid with periodic boundary condition. Let $`L`$, the linear dimension of the system (i.e. $`N=L^2`$), be even. For a non-zero parameter value we (by uniform randomness) delete $`r_1L^2`$ ‘diagonal’ edges creating frustration as illustrated in Fig. 2(b). To be more precise, if we index the vertices as $`(i_x,i_y)`$, $`1i_x,i_yL`$; then the edges are $`[(i_x,i_y),(i_x+1,i_y)]`$ and $`[(i_x,i_y),(i_x,i_y+1)]`$ (giving the square grid) plus $`r_1L^2`$ edges of the form $`[(i_x,i_y+1),(i_x+1,i_y)]`$ chosen by uniform randomness (addition is modulo $`L`$). This model has a high degree of short circuits. The extremes $`r_2=0`$ and $`r_1=1`$ represent two generic lattice types. The symmetries of the regular networks simplify the calculations of e.g. limit properties for the bipartivity measures. If $`r_2=1`$ the system is bipartite (note that $`L`$ has to be even for this to hold) so $`b_{1,2}=1`$. When $`r_2=0`$ we have $`b_1=b_2=2/3`$: For the lower limit of the $`b_1`$ quantity, see Ref. WAHO . For the lower limit $`b_2`$ we note that $`\mathrm{\Sigma }(C_3)=6N`$ (since each vertex can be associated with two triangles). This gives $`\widehat{n}=3`$ and $`\nu =2`$ for all edges. Now it is enough to mark $`N`$ edges (e.g. all $`[(i_x,i_y),(i_x+1,i_y)]`$ edges). In this case we note that each edge will have $`\nu =2`$ when it is marked, which means that the marking sequence is optimal and that the number of iteration cannot be less with another choice of edges to mark. So $`b_2=1N/3N=2/3`$. The major disadvantage with Model 2 is that the average degree is a function of $`r_2`$ ($`M=(3r_2)L^2`$). This change in the average degree can make it harder to separate effects of the shift in bipartivity from the shift in average degree.

In both model 1 and (even more) model 2 triangles will dominate the set of odd circuits. To test networks with predominantly longer circuits we construct a Model 3 as follows (see Fig. 2(c)): We make two circulants of size $`N/2`$ with the vertices $`\{v_1^i,\mathrm{},v_{N/2}^i\}`$ and edges $`\{(v_1^i,v_2^i),\mathrm{},(v_{N/21}^i,v_{N/2}^i),(v_{N/2}^i,v_1^i)\}`$, $`i\{1,2\}`$. Then we add $`M_{\mathrm{trans}}`$ transverse edges between the circulants. $`M_{\mathrm{trans}}/2`$ of these edges are placed out separated by equal distance $`N/M_{\mathrm{trans}}`$ separating the double circulants into $`M_{\mathrm{trans}}/2`$ ‘sectors.’ Then we fill up each sector with another transverse edge: With probability $`r_3`$ we add an $`(v_i^1,v_i^2)`$ edge (such that $`(v_i^1,v_i^2)`$ is none of the previously added transverse edges), otherwise we add a $`(v_i^1,v_i^2+1)`$ edge (addition modulo $`N/2`$). We note, to a first approximation, that if $`r_3=0`$ marking (in the process of calculating $`b_2`$) one edge between every transverse edge on one of the circulants is needed to mark the shortest odd circuits. This will make $`b_2O(1M_{\mathrm{trans}}/N)`$.

### III.2 Real-world networks

Physicists’ networks studies has, in the spirit of statistical mechanics, emphasized the properties remaining when the system grows beyond any limit. Bipartivity, as discussed above, is well defined for all system sizes. Still it is a quantity that can potentially suffer from finite-size effects (from the fact that not all real neighbors of all actors in a empirically constructed social network are a part of the graph) and is therefore preferably measured for large networks. Now the problem is to find data for large-scale real-world networks of social interaction. In general two methods has been successful for this purpose—one either uses professional collaborations of some sort or data from interaction over the Internet (either in Internet communities HEL ; smith , or through email exchange ebel .

#### III.2.1 Professional collaboration networks

In the professional collaboration networks we study the vertices are professionals of some field—networks of scientists and company directors are considered in this papers, the movie-actor network is another frequently studied example; the edges represent that two actors has been involved in the same professional collaboration. This is some-times referred to as a “one-mode” representation of an affiliation network (as opposed to the bipartite two-mode representation discussed in Sect. I).

Professional collaboration networks are no doubt interesting in their own right as accounts for the interaction dynamics of the respective fields. Assuming that the formation of professional ties follow similar principles as general human interaction, we can use professional collaboration networks to draw conclusions about the structure of more general social networks. However, at one point (at least) professional collaboration differs from general social interactions: A collaboration tie does not necessarily imply a strong personal acquaintance, but in these networks each collaboration constitutes a fully connected cluster. This leads to higher fraction of short circuits than, say, a friendship network.

One of the professional collaboration network we use is of scientists who has uploaded manuscripts to the preprint repository arxiv.org. Two scientists are linked if their name (identified by surname and initials) appear together on at least one preprint. A detailed description of this network can be found in Ref. newman3 . In the other professional collaboration network the vertices represent company directors from the Fortune top 1000 list of companies in USA the year 2001. An edge (collaboration) in this network means that two directors are sitting in board of the same company. A detailed description of this network can be found in Ref davis . Sizes of the networks can be seen in Table 1.

#### III.2.2 Online interaction networks

In online interaction networks, the vertices are users of Internet communities and an arc (A,B) is added if A contacts B, or if A adds B to his/her list of friends smith ; HEL . Another kind of online interaction networks are email networks ebel , where an arc can be assigned if an email is sent, or if a person adds another to his/her address book. Just as for professional collaboration networks, one can argue that online interaction networks are representative as general social networks. One can assume that new contacts are formed through preference-matching searches to a larger extent, and introduction by mutual friends to a lesser extent, than in general friendship networks. Since the introduction of mutual friends to each other is believed to be the major cause of high clustering (large density of triangles, or, large transitivity) newman1  one can expect a lower clustering in networks of online interaction (still the clustering in these network seems to be finite in the $`N\mathrm{}`$ limit HEL ).

The specific online interaction networks we consider are constructed from the Internet communities nioki.com and pussokram.com. The nioki.com data is described in Ref. smith . In this data an arc (A,B) means that B is listed as a friend by A, which allows A to see if B is online and send instant messages to B. In the pussokram.com data the arcs correspond to communication between the users. There are four different types of communication in this specific network (all described in detail in Ref. HEL ). We use the networks obtained from two types of interaction (‘messages’—like ordinary emails within the community, and “guest book”—where one user contacts another by writing in his/her guest book), and the network of any of the four types. Network sizes can be found at Table 1.

Another large difference between the pussokram.com and nioki.com data is that the former community has a very pronounced romantic profile, encouraging flirts and romantic correspondence. nioki.com has also a search engine to “trouve l’amour” (find love), but that is all.

Apart from the two Internet communities, we study another type of online interaction network based on the flow email. For this network all in- and out-going email traffic to a server was logged for around three months ebel . The server handles undergraduate students’ email accounts at Kiel University, Germany. Thus there are two categories of vertices—internal vertices, whose activity is accurately mapped; and external vertices, that only have edges leading to internal vertices. In this study we restrict ourselves to the network of internal-internal contacts. The reason we do not include external contacts is that we would miss the (probably many) circuits containing external-external edges which would bias the bipartivity.

#### III.2.3 Network from interview and field survey

Apart from the above networks, all obtained from databases, we also measure the bipartivity of two networks obtained from interview and field surveys. The first data set is gathered by observations of interaction between members of a university karate club zach . We also study the network of acquaintance ties in a prison prison . The outgoing arcs from A corresponds to prisoners listed by A in response to the question: “What fellows on the tier are you closest friends with?” Due to their acquisition methods these kind of real-world networks has to be rather small. This can, as mentioned, result in finite size effects. On the other hand they, most likely, more truly reflect the structure of real acquaintance networks.

## IV Results

In this section we present the results of the test networks and the measurement for the real-world social networks.

### IV.1 Test networks

As expected, both $`b_1`$ and $`b_2`$ are monotonously increasing as functions of the $`r_1`$, $`r_2`$ and $`r_3`$ parameters of (almost note:not\_really  all our test network (see Fig. 3)). This is encouraging and suggests that both $`b_1`$ and $`b_2`$ are quite relevant measures of bipartivity.

The Model 1 measurements shown in Fig. 3(a) are made with the model parameters $`N=2\stackrel{~}{N}=100`$ and $`M=800`$. We have checked many other sizes too, but all have the characteristic appearance of Fig. 3(a)—a linear increase of $`b_1`$ and $`b_2`$ for larger $`r_1`$ and an flatter slope for $`r_1`$ close to zero. This shape is expected from the discussion in Sect. I—in networks where a heterophilous preference is the only structure-inducing force, only the strong preference limit gives a strong measurable effect: Close to the ER limit $`r_10`$, the original two partitions will not be identified correctly, only when the different partition (to a large extent) have different sign the bipartivity will be proportional to the strength of the heterophilous preference.

As seen in Fig. 3(b) Model 2 shows an almost linear functional form of $`b_{1,2}(r_2)`$. In this case, triangles dominate the odd circuits even at small values of $`r_2`$. Tuning $`r_2`$ will give a proportional increase of the number of triangles. Thus a linear $`r_2`$ dependence of $`b_2`$ would be expected.

Also Model 3 has linear $`b_{1,2}`$ vs. $`r_3`$ curves. The model parameters used are $`N=100`$ and $`M_{\mathrm{trans}}=10`$. As mentioned in Section III.1, we expect $`b_2M_{\mathrm{trans}}/N`$ for $`r_3=0`$, which is confirmed in Fig. 3(c).

The measurements for both $`b_1`$ and $`b_2`$ are averaged over 100 network realizations. The XMC scheme for the $`b_1`$ quantity is ran at 24 temperatures in parallel, between temperatures $`0.01`$ and $`2`$. Other network parameters are $`t_{\mathrm{avg}}=4\times 10^5`$, $`t_{\mathrm{measure}}=4`$, $`t_{\mathrm{quench}}=20`$ and $`t_{\mathrm{exch}}=1000`$. These are more modest parameter values than we will use for the real-world networks, but the test networks are also much smaller, and since the distribution of $`b_1`$ and $`b_2`$ are (likely) symmetric, the network average helps to reduce the error.

### IV.2 Real-world social networks

Now we turn to the result for the bipartivity measures of real-world networks. The values are presented in Table 1. For comparison we also give values for the clustering coefficient (density of triangles) $`C`$ and the density of squares $`D`$ in both directed and undirected versions note:cd . Undirected networks are constructed by taking the reflexive closure. At first glance at the table we arrive at the pleasing conclusion that the bipartivity for the pussokram.com networks is very high (as expected from a network of romantic interaction of mostly heterosexuals). But disappointingly, the bipartivity measures show similarly high values for the nioki.com and email networks. This can be explained by the fact that nioki.com, just like the pussokram.com, data has very low $`C`$ and $`D`$ values, and presumably very few circuits at all. Now branches (subgraphs without circuits that can be isolated by cutting one edge) does not give a positive contribution to either $`b_1`$ or $`b_2`$, no matter of the gender of the agents. The email network do have a high clustering, but still rather high bipartivity. The reason is that the email network is rather heavily fragmented and contains many isolated subnetworks of two vertices and one edge, and three vertices and two edges. Such subnetworks does not affect the clustering coefficient but tends to decrease the bipartivity measures note:improvement .

The collaboration networks consist of a number of fully connected clusters (corresponding to a specific collaboration) that are interconnected. It is thus natural that we see low bipartivity and a high density of short circuits. The lower bipartivity values for the company director network can be explained by smaller average size of such fully-connected clusters: The average number of vertices per collaboration is $`9.5`$ for the corporate director network and $`2.5`$ for the scientific collaboration data davis ; newman3 .

The two small networks constructed from field surveys (the “karate club” and “prison” network of Table 1, discussed in Section III.2.3) show mid-range bipartivities and relative high values of $`C`$ and $`D`$. From the above discussion we can expect that the bipartivity of large, real, acquaintance networks is somewhere between those of the collaboration networks and the Internet community networks (because they probably have higher clustering than Internet community networks, and lower number of fully connected clusters than the collaboration networks). Encouraging enough, this is exactly what we see in Table 1. Of course, the very small systems sizes might affect the results, but that the bipartivity measures of real-world acquaintance measures would be close to either the upper or lower limits seems hard to believe.

We conclude this section by a note on the parameters for the XMC optimization. The measurement of $`b_1`$ for all real-world network (except the nioki.com data where we study the convergence more carefully) are done just once with the following simulation parameters, $`N_T=24`$ (with temperatures from $`0.002`$ to $`5`$) $`t_{\mathrm{avg}}=1\times 10^7`$, $`t_{\mathrm{measure}}=16`$, $`t_{\mathrm{quench}}=40`$ and $`t_{\mathrm{exch}}=2\times 10^4`$.

## V Summary and discussion

This paper concerns the quantification of the network structure ‘bipartivity’—how close to bipartite a given graph is. We propose two measures for this quantity. One quantity $`b_1`$ based on the optimal two-coloring of the network—or, equivalently, the ground state of the antiferromagnetic Ising model on the network. The exact value of this quantity (that has been used in different roles elsewhere) is NP-complete and thus in general not feasible to calculate exactly. Instead we seek an approximate solution by a simulated annealing approach. The simulated annealing is based on the exchange Monte Carlo scheme. We argue that this unorthodox minimization method helps us avoid local minima of the energy landscape of the antiferromagnetic Ising model. Furthermore we develop a measure $`b_2`$ based on the count of odd circuits that, for almost all networks, is calculable in polynomial time.

We propose three different random graph test models where one can interpolate between arguably non-bipartite and bipartite graphs by tuning a control parameter. Both our bipartivity measures are shown to increase monotonically with tuning the control parameters towards the bipartite extreme. From this we conclude that the bipartivity measures really quantify the notion of bipartivity.

By considering example networks we infer that bipartivity is a structure that cannot be measured by currently popular structural measures, such as the clustering coefficient. At the same time any sensible quantification of bipartivity probably has to have a positive correlation with the clustering coefficient for most networks (with exceptions for exotic cases like Fig. 1(a))—so, in that case bipartivity and clustering is not independent.

We measure $`b_1`$ and $`b_2`$ of a number of real-world networks, constructed from online interaction, professional collaborations, and field surveys. As expected, we see high bipartivity values for data from the Internet community pussokram.com, where romantic contacts are encouraged, and hence a high degree of heterophilous interaction expected. We also see the expected low bipartivity values for the professional collaboration and empirical acquaintance networks we study. Disappointingly we cannot use our bipartivity measures to distinguish between the networks driven by romantic or friendship (or professional) contacts. To do this other structures and the network sizes has to be taken into account, in a more elaborate analysis (that is out of the scope of this study).

So far our examples of networks with high bipartivity has been romantic networks and networks of sexual contacts. Network-based studies of sexually transmitted diseases lea  is a potentially interesting area for bipartivity measures, as the transmission rates for homosexual and heterosexual contacts differ anma . Apart from romantic and sexual networks, there are other areas where the bipartivity measure may prove useful: One can consider a trade network where some agents are more or less pronounced sellers and others are primarily buyers (cf. Ref. white ), such networks would not have a neutral bipartivity. Another application is for the ‘genealogical’ network of a disease outbreak: Some contagious diseases have a relatively stable duration between when an individual is infected and when he or she becomes infectious. Epidemics of these types of diseases can therefore roughly be divided into different generations of infected individuals anma . A network consisting of possible edges of infections, for an outbreak of this type of disease, should therefore have very few odd-length circuits. The reason is that the infection is only transmitted between succeeding generations, which generates only circuits of even length (in the reflexive closure of the network). When reconstructing the paths this kind of disease has taken in a population, a minimization of the bipartivity measures can be a method for excluding redundant infectious edges.

We conclude by an analogy to linear algebra—we have identified a new dimension (structure) and proposed base vectors (measures), that unfortunately are not orthogonal to the other dimensions.

## Acknowledgements

We would like to thank Niklas Angemyr, Stefan Bornholdt, Gerald Davis, Holger Ebel, Michael Lokner, Stefan Praszalowicz, and Christian Wollter for help with data acquisition; and Johan Giesecke, James Moody, Mats Nylén, and Pontus Svenson, for comments and suggestions. P.H. was partly supported by the Swedish Research Council through contract no. 2002-4135. F.L. was supported by the National Institute of Public Health. C.R.E. was supported by the Bank of Sweden Tercentenary Foundation. B.J.K. was supported by the Korea Science and Engineering Foundation through Grant No. R14-2002-062-01000-0.

## Appendix A The lower bound of the measure $`b_2`$

In this Appendix we argue that, in the $`N\mathrm{}`$ limit, the lower bound for $`b_2`$ is $`1/2`$ (just like $`b_1`$). First we conjecture that the minimal value for $`b_2`$, just as for $`b_1`$, is attained for complete graphs. (This will be further motivated below.)

To assess $`b_2`$ for complete graphs, we note that note:sigsum

$`\mathrm{\Sigma }(C_n)`$ $`=`$ $`{\displaystyle \underset{\mathrm{odd}\mathrm{\hspace{0.33em}3}in}{}}{\displaystyle \frac{N!}{2(Ni)!}}`$ (11a)
$`\mathrm{\Sigma }(C_3)`$ $`=`$ $`{\displaystyle \frac{N(N1)(N2)}{2}}`$ (11b)
$``$ $`{\displaystyle \frac{N(N1)}{2}}=M,`$
so $`\widehat{n}=3`$ which results in that $`\nu =N2`$ for each edge.

Now we apply the marking procedure of Sec. II.2.1. Marking an edge $`(u,v)`$ makes $`\nu (u,v)=\nu (v,u)=0`$. Furthermore, every edge $`(u,w)`$ and $`(v,w)`$ ($`wu,v`$) will be decreased by one since the triangle $`\{u,v,w\}`$ now contains a marked edge. The discussion will be simplified by considering a matrix representation of $`\nu (u,v)`$. Marking $`(u,v)`$ sets $`\nu (u,v)=\nu (v,u)=0`$ and decreases the $`u`$’th and $`v`$’th columns, and $`u`$’th and $`v`$’th rows by one (an example is given in Fig 4(a)). Marking another edge $`(u^{},v^{})`$ ($`u^{}`$ and $`v^{}`$ are different from both $`u`$ and $`v`$, otherwise $`\nu (u^{},v^{})`$ would not be maximal) will have the same effect as marking the first. For positions like $`(u,v^{})`$ the original $`\nu `$ are decreased by 2 (see Fig. 4(b)), since it has lost the two passing triangles $`\{u,u^{},v^{}\}`$. and $`\{v^{},u,v\}`$. Continuing this process we see that it takes $`N/2+O(1)`$ markings for $`\nu `$ of each edge to be decreased by two units, and thus $`m^{}=N^2/4+O(N)`$ markings to make $`\nu =0`$ for all edges. This gives $`b_2=1/2`$ in the $`N\mathrm{}`$ limit. Since the appropriateness of $`b_2`$ as a bipartivity measure is not really dependent on the limit values, we will not give a rigorous proof that the correction is of a lower order for all levels of the marking procedure (one level is the $`N/2+O(1)`$ edges needed to be marked for $`\nu `$ to be decreased by at least two units for each edge).

Now we argue that the $`b_2`$ takes its minimal value for complete graphs. First we note that the number of circuits of length $`n`$ per edge, for any $`n`$, is largest in a complete graph review2 . So if we set $`\widehat{n}`$ arbitrarily and discard circuits of length $`n`$ in the calculation of $`\nu (v)`$, the fully connected graph would give the highest $`m^{}`$ value and thus the lowest bipartivity measure. The strongest candidate for a lower bipartivity measure than that of a fully connected graph would thus be a graph such that the $`\mathrm{\Sigma }(C_n)<3M`$ and $`\mathrm{\Sigma }(C_{n+2})`$ is as big as possible for some $`n`$. But the number edges needed to be removed from a fully connected graph for $`\mathrm{\Sigma }(C_n)<3M`$ to hold, not only reduces the contribution to $`\nu `$ from circuits of length $`n`$ but also from circuits of length $`n+2`$ to a similar extent. If one performs the approximate marking procedure outlined above for circuits of length five one starts from $`\nu =(N2)(N3)(N4)`$ and it takes $`N/2+O(1)`$ markings to decrease every $`\nu `$ with at least $`2N^2`$. This means that the number of edges needed to be marked to make $`\nu =0`$ for every edge is the same if circuits of length five is considered. It also means that a graph as outlined above (with $`\mathrm{\Sigma }(C_n)<3M`$ and $`\mathrm{\Sigma }(C_{n+2})`$ is as big as possible) probably do not have a lower $`b_2`$ than a complete graph.

To epitomize, the $`b_2`$ measure lies in the interval $`[1/2,1]`$ in the $`N\mathrm{}`$ limit. The finite size corrections to $`b_2`$ for fully connected graphs, however, turns out to make $`b_2`$ slightly less than $`1/2`$.

## Appendix B Convergence of the simulated annealing

To analyze the convergence of the simulated annealing scheme we run ten independent calculations of the $`b_1`$ quantity (with the same parameter values as in Sect. IV.2). The individual time evolutions of $`b_1`$ (at the lowest temperature $`T=0.002`$) for the different runs are shown in Fig. 5. We note that already after the first quench $`b_1`$ is only $`3\%`$ away from the value at the end of the run, and after 50 time steps $`b_1`$ is $`0.5\%`$ of the value after $`1\times 10^7`$ time steps. We note that there is no way of constructing a statistically valid confidence interval for the true $`b_1`$ value since an arbitrary complex energy landscape could have a global minimum with a basin of attraction of measure zero. There are however indications that this is seldom a major problem, at least not for the bisection problem jerrum .

An interesting observation from Fig. 5 is the step-like structure. This is a result of the exchange trials: After $`t100`$ the local minimum has been found, but at the temperature in question the system is in principle stuck in a confined part of the configuration space, and cannot enter lower lying energy valleys. In the time scale $`t=10^5`$ there is another jump in the $`b_1`$ value. This is related to that other replicas from other parts of the configuration space reaches the lowest level. At around $`t=10^6`$ the current highest $`b_1`$ values (lowest energy) reaches another plateau. At this time, each replica should have covered the whole temperature range several times. This second plateau gives two encouraging implications: Firstly, that the correct value of $`b_1`$ probably is not very far off the measured value. Secondly, that the exchange steps really are helpful. If one wants to run this algorithm more efficiently the $`t_{\mathrm{exch}}`$ we use is far too large (but beneficial for separating the time scales in the discussion above). Ideally $`t_{\mathrm{exch}}`$ should probably be chosen to be of the same order as the first jump (from the regular Monte Carlo steps)—in the nioki.com network (displayed in Fig. 5) this would be $`t100`$.


text:
# Phase-field approach for faceted solidification

## I Introduction

Over the last decade, the phase-field approach  has been developed extensively to model the solidification of both pure materials and alloys . Most of the work to date has focused on the case where the excess free-energy of the solid-liquid interface

$$\gamma \gamma _0f(\theta )$$
(1)
is a smooth function of the angle $`\theta `$ between the direction normal to the interface and some fixed crystalline axis. In particular, the simple form

$$f(\theta )=1+ϵ\mathrm{cos}4\theta +\mathrm{},$$
(2)
appropriate for a weakly anisotropic material with an underlying cubic symmetry has been widely used in studies of dendritic solidification . For a smooth $`\gamma `$-plot, the value of the diffusion field (dimensionless undercooling or supersaturation) at the interface is given by the standard Gibbs-Thomson condition

$$u=d_0\left[f+\frac{d^2f}{d\theta ^2}\right]\kappa ,$$
(3)
where $`d_0`$ is a microscopic capillary length (thermal or chemical) that is proportional to $`\gamma _0`$ and $`\kappa `$ is the interface curvature. Dendritic growth in pure materials  and alloys  has been modeled quantitatively with a thin-interface limit  of the phase-field model that yields the boundary condition (3) and overcomes the stringent computational constraint associated with a finite interface thickness. These simulation studies modeled the common situation of a weak anisotropy where the interface stiffness $`\gamma +d^2\gamma /d\theta ^2`$, and hence the square bracket in Eq. (3), is everywhere positive. More recently, the phase-field approach has also been successfully extended to model strongly anisotropic equilibrium crystal shapes where the stiffness becomes negative for some range of $`\theta `$ . In this case, which occurs for example when $`ϵ>1/15`$ for the form (2) of $`f(\theta )`$, the slope of the interface has discontinuities that must be present to eliminate thermodynamically unstable, and hence forbidden, orientations for which the stiffness is negative.

In this paper, we extend the phase-field approach to model a wide class of materials that form facets for a discrete set of orientations. The interfacial energy is generally non-analytic for orientations close to a facet. This non-analyticity is reflected in the presence of cusps in the $`\gamma `$-plot that are of the form

$$f(\theta )1+\delta |\theta \theta _c|+\mathrm{}\mathrm{for}|\theta \theta _c|1,$$
(4)
where $`\theta =\theta _c`$ is the orientation of a given facet. In the simplest model where vicinal surfaces are assumed to consist of straight ledges of height $`a`$ separated by terraces of width $`la/|\theta \theta _c|`$, the cusp amplitude $`\delta =\gamma _L/(a\gamma _0)`$ where $`\gamma _L`$ is the ledge energy per unit length.

To a good first approximation, the tendency of a material to facet is described by Jackson’s well-known $`\alpha `$-factor  that is the product of a dimensionless crystallographic factor that depends on the crystal structure and the orientation of the interface and the ratio of the latent heat per mole to the rare gas constant. Jackson’s theory predicts that materials with $`\alpha 2`$ (such as metals with low entropy of melting) will have rough interfaces, while those with $`\alpha 2`$ (such as non-metals and compounds) will form facets, consistent with experimental observations across a wide range of materials .

The growth of faceted dendrites has been studied both experimentally  and theoretically in the context of a sharp-interface approach . Experimentally, facets are seen to appear near the tip of the needle crystal that nevertheless retains a parabolic shape on the larger scale of the whole crystal that includes the trailing rough parts. From a theoretical standpoint, a main difficulty to model faceted growth is that the standard Gibbs-Thomson relation can only be applied to rough parts of the interface. On a facet, this relation takes the form of an integral condition that is obtained by making the substitution $`\kappa =d\theta /ds`$, where $`s`$ denotes the arclength along the interface, and integrating both sides of Eq. (3) from one extremity ($`s_{}`$) of the facet to the other ($`s_+`$). This integration yields the condition

$$_s_{}^{s^+}u𝑑s=d_0\left(\frac{df}{d\theta }|_+\frac{df}{d\theta }|_{}\right)=2d_0\delta ,$$
(5)
where $`df/d\theta |_\pm `$ denotes the limiting values of $`df/d\theta `$ on each side of the cusp ($`\theta \theta _c0^\pm `$).

For an equilibrium crystal with a constant undercooling or supersaturation $`(u<0)`$ along the interface, the above condition implies that the length $`\mathrm{\Lambda }`$ of the facet is simply proportional to the amplitude of the cusp

$$\mathrm{\Lambda }=2\frac{d_0}{u}\delta .$$
(6)
For growth outside of equilibrium, the determination of the facet length and the shape of the rough parts requires in general a self-consistent solution of the free-boundary problem defined by the appropriate diffusion equations for $`u`$ in solid and liquid, the standard Stefan condition of heat or mass conservation on the interface, and the local equilibrium conditions (3) and (5) on rough parts and facets, respectively. We defer a discussion of kinetic effects to the conclusion section of this paper.

An analytical solution to this problem was obtained for steady-state growth by Adda-Bedia and Hakim  in the small Péclet number limit neglecting capillary effects on the rough parts. They concluded from this analysis that it is not possible to require both tangential matching of the rough and faceted parts of the interface and the equilibrium condition ($`u0`$) on the front and trailing rough parts. Furthermore, they proposed an approximate solution to the full problem that includes capillary effects on the rough parts in the large $`\delta `$ limit. In this solution, the facets are assumed to match tangentially to a small quasi-circular tip that is significantly more undercooled than the rough trailing parts. A numerical solution of the steady-state growth problem was later obtained for arbitrary $`\delta `$ by Adda-Bedia and Ben Amar  using a boundary integral method. This calculation gave results that are consistent with the the approximate solution of Adda-Bedia and Hakim for large $`\delta `$.

At present, a phase-field approach for faceted growth would be highly desirable to study the stability of needle crystal solutions as well as to explore the full dynamical range of morphological evolution when the interface shape is non-stationary. Treating separately the rough and faceted parts of the interface appears to be difficult within a phase-field formulation. To circumvent this difficulty, we follow a procedure that is similar in spirit, but different in details, to the one introduced by Adda-Bedia and Ben Amar  in a sharp-interface context. The basic idea is to use a regularized form of the $`\gamma `$-plot that can approximate arbitrarily closely the $`\gamma `$-plot with sharp cusps. The discontinuity of $`df/d\theta `$ at a cusp can be viewed mathematically as the stiffness behaving as a sharply peaked Dirac delta function at $`\theta =\theta _c`$. Therefore, different regularizations of this discontinuity can be viewed as different regularizations of this delta function peak that must satisfy the cusp condition

$$\underset{\theta _00}{lim}_{\theta _c\theta _0}^{\theta _c+\theta _0}𝑑\theta \left[f+\frac{d^2f}{d\theta ^2}\right]=2\delta ,$$
(7)
imposed by the form (4) of $`f`$ near a cusp. Adda-Bedia and Ben Amar have introduced regularizations in which all derivatives of the stiffness and hence $`f`$ are continuous. In a phase-field context, we have found more convenient to use a regularization where only $`f`$ and its first derivative are continuous and the dimensionless stiffness $`f+d^2f/d\theta ^2`$ is a step function regularization of the delta function, $`f+d^2f/d\theta ^2\delta /\theta _0`$ for $`\theta _c\theta _0\theta \theta _c+\theta _0`$. In practice, computations for small $`\theta _0`$ turn out to be more costly because the larger stiffness ($`\delta /\theta _0`$) requires a smaller time step. As we shall see, however, an accurate extrapolation to the $`\theta _00`$ limit is nevertheless possible.

The phase-field model and the details of this smoothing procedure are described in the next section. The numerical implementation of the equations is then presented in section III. The convergence of our numerical approach is demonstrated in section IV by a detailed comparison of diffuse-interface and sharp-interface equilibrium crystal shapes. Results concerning the steady-state growth of faceted needle crystals in a pure undercooled melt are then presented in section V. These results are then compared in section VI to the predictions of an approximate sharp-interface theory that assumes simple forms for the front and trailing rough parts as an extension of the solution proposed in Ref. . Concluding remarks and future prospects are then discussed in section VI.

## II Phase-field Model

We consider the growth of a crystal from a pure undercooled melt with the thermal diffusivity being the same in the solid and liquid. We follow the phase-field approach developed by Karma and Rappel  that has proven successful to model efficiently non-faceted dendritic growth. We briefly recall the basic features of the phase-field model and of the thin-interface limit used to relate the phase-field and sharp-interface models.

The dimensionless temperature field is defined as $`u=(TT_M)/(L/c_p)`$, where $`T`$ is the temperature, $`T_M`$ the melting temperature, $`L`$ the latent heat per unit volume, and $`c_p`$ the specific heat at constant pressure per unit volume. The phase field $`\psi `$, is taken to be equal to $`+1`$ ($`1`$) in the solid (liquid) phase, and varies continuously across the diffuse interface. In two dimensions, the time evolution of both fields is governed by the equations

$$_tu=D^2u+\frac{1}{2}_t\psi ,$$
(8)
and

$`\tau (\theta )_t\psi `$ $`=`$ $`\left[\psi \lambda u(1\psi ^2)\right](1\psi ^2)`$ (11)
$`+\stackrel{}{}\left[W(\theta )^2\stackrel{}{}\psi \right]_x\left[W(\theta )W^{}(\theta )_y\psi \right]`$
$`+_y\left[W(\theta )W^{}(\theta )_x\psi \right],`$
where $`D`$ is the thermal diffusivity, $`\theta `$ is the angle between the normal to the interface and the $`x`$ axis, $`\lambda `$ is a dimensionless coupling constant between the temperature and phase fields, and

$$W(\theta )W_0f(\theta ),$$
(12)
is the diffuse interface thickness. In the thin-interface limit where $`W_0`$ is much smaller than the mesoscale of the growing crystal , the above phase-field equations reduce to the standard sharp-interface model of diffusion-limited growth with the velocity-dependent form of the Gibbs-Thomson condition

$$u=d_0\left[f(\theta )+\frac{d^2f(\theta )}{d\theta ^2}\right]\kappa \beta (\theta )v_n,$$
(13)
where

$$d_0=\frac{a_1W_0}{\lambda },$$
(14)
and

$$\beta (\theta )=\frac{a_1}{\lambda }\frac{\tau (\theta )}{W(\theta )}\left[1a_2\lambda \frac{W(\theta )^2}{D\tau (\theta )}\right].$$
(15)
In addition, $`a_1`$ and $`a_2`$ are constants that depend generally on the choice of the double-well potential and other functions in the phase-field model. For the present choices that are the same as in Refs. , $`a_1=0.8839\mathrm{}`$ and $`a_2=0.6267\mathrm{}`$.

To model faceted growth, we use the simplest form of the $`\gamma `$-plot

$$f(\theta )=1+\delta \left(|\mathrm{sin}\theta |+|\mathrm{cos}\theta |\right).$$
(16)
This form is directly adapted from the broken-bond model which describes satisfactorily faceted solid-gas interfaces . The same form was used in previous sharp-interface calculations .

Furthermore, in order to model the limit with local equilibrium at the interface, i.e vanishing interface kinetics $`\beta (\theta )=0`$, we impose

$$\tau (\theta )=\tau _0f(\theta )^2,$$
(17)
and

$$\lambda =\frac{1}{a_2}\frac{D\tau _0}{W_0^2},$$
(18)
which makes $`\beta (\theta )`$ vanish for all values of $`\theta `$.

For the above choice of $`\gamma `$-plot, the dimensionless stiffness defined by

$$S(\theta )f(\theta )+\frac{d^2f(\theta )}{d\theta ^2}$$
(19)
is constant except at the cusps, $`\theta =\theta _{c,n}=n\pi /2`$ (with $`n`$ an integer), where it diverges. The discontinuities in the first derivative of the interfacial energy need to be regularized, since $`f^{}(\theta )`$ is not defined for $`\theta =\theta _{c,n}`$. A simple strategy to circumvent this problem is to smooth out the cusps by replacing $`f(\theta )`$ with a smooth function like a sine in a small range of $`\theta `$ values around the cusps (Fig. 1). For the sake of clarity, we restrict our discussion to the first quadrant, $`\theta [0,\frac{\pi }{2}]`$, since the problem has a natural four-fold symmetry. The smoothed anisotropy function $`f_s(\theta )`$ reads then:

$$f_s(\theta )=\{\begin{array}{cc}1+\delta \left(\mathrm{sin}\theta +\mathrm{cos}\theta \right)\hfill & \mathrm{if}\theta _0<\theta <\frac{\pi }{2}\theta _0,\hfill \\ BA\mathrm{cos}\theta \hfill & \mathrm{if}\theta \theta _0,\hfill \\ BA\mathrm{sin}\theta \hfill & \mathrm{if}\theta \frac{\pi }{2}\theta _0,\hfill \end{array}$$
(20)
and its derivative

$$f_s^{}(\theta )=\{\begin{array}{cc}\delta \left(\mathrm{cos}\theta \mathrm{sin}\theta \right)\hfill & \mathrm{if}\theta _0<\theta <\frac{\pi }{2}\theta _0,\hfill \\ A\mathrm{sin}\theta \hfill & \mathrm{if}\theta \theta _0,\hfill \\ A\mathrm{cos}\theta \hfill & \mathrm{if}\theta \frac{\pi }{2}\theta _0.\hfill \end{array}$$
(21)
The two constants $`A`$ and $`B`$ are obtained by expressing the continuity of $`f_s`$ and $`f_s^{}`$ at $`\theta =\theta _0`$,

$$A=\delta (\mathrm{cot}\theta _01),$$
(22)
and

$$B=1+\delta /\mathrm{sin}\theta _0.$$
(23)
Now $`f_s^{}(0)=f_s^{}(\pi /2)=0`$ and $`W(\theta )`$ is no more singular in the cusp directions. Note, however, that $`f_s^{\prime \prime }`$ is not continuous at $`\theta =\theta _0`$. Consequently, the dimensionless stiffness is a step function of $`\theta `$ (see Fig. 1),

$$S(\theta )=\{\begin{array}{cc}1\hfill & \mathrm{if}\theta _0<\theta <\frac{\pi }{2}\theta _0,\hfill \\ 1+\delta /\mathrm{sin}\theta _0\hfill & \mathrm{if}\theta \theta _0\mathrm{or}\theta \frac{\pi }{2}\theta _0.\hfill \end{array}$$
(24)
## III Numerical implementation

We now briefly describe how Eqs. (8, 11) are discretized in our code. The interface is represented by the $`\psi =0`$ contour, so that $`\stackrel{}{}\psi `$ is by construction colinear to the unit vector along the normal to the interface, n. If $`|\stackrel{}{}\psi |0`$, the two components of this vector are given by

$$n_x=\mathrm{cos}\theta =_x\psi /|\stackrel{}{}\psi |,$$
(25)
and

$$n_y=\mathrm{sin}\theta =_y\psi /|\stackrel{}{}\psi |,$$
(26)
and $`f_s`$ and $`f_s^{}`$ are computed according to Eqs (20, 21). Conversely, when $`|\stackrel{}{}\psi |=0`$, one sets

$$f_s(\theta )=1,$$
(27)
and

$$f_s^{}(\theta )=0.$$
(28)
Numerical simulations of Eqs. (8, 11) are performed by implementing a finite-difference scheme on a regular square mesh with mesh size $`h=0.4W_0`$. The domain considered is the quadrant $`0xNh`$ and $`0yNh`$, with $`N`$ an integer. A standard first-order in time Euler scheme with a time step $`\mathrm{\Delta }t`$ and a second-order in space discretization of the spatial derivatives are used.

Let $`u_{i,j}`$ and $`\psi _{i,j}`$ respectively denote the discretized reduced temperature and phase field at point $`(x=ih,y=jh)`$. In a first step, one computes $`^2u`$, $`^2\psi `$, $`_x\psi `$, $`_y\psi `$, $`W(\theta )W^{}(\theta )`$ and $`W(\theta )`$ for each point of the extended domain, $`i`$ and $`j[1,N+1]`$, and the results are stored in six intermediate arrays. Reflecting conditions are imposed on $`u`$ and $`\psi `$ at all the domain boundaries ($`i`$ or $`j=2,1,N+1,N+2`$), and using centered differencing approximations, we have

$$^2u_{i,j}=\left[u_{i+1,j}+u_{i1,j}+u_{i,j+1}+u_{i,j1}4u_{i,j}\right]/h^2,$$
(29)
$$^2\psi _{i,j}=\left[\psi _{i+1,j}+\psi _{i1,j}+\psi _{i,j+1}+\psi _{i,j1}4\psi _{i,j}\right]/h^2,$$
(30)
$$_x\psi _{i,j}=\left[\psi _{i+1,j}\psi _{i1,j}\right]/(2h),$$
(31)
$$_y\psi _{i,j}=\left[\psi _{i,j+1}\psi _{i,j1}\right]/(2h).$$
(32)
The two remaining arrays, $`(WW^{})_{i,j}`$ and $`W_{i,j}`$, are given by Eq. (12) into which $`f`$ is replaced with $`f_s`$.

In the second step, one solves Eqs. (8) and (11) on the inner points of the domain, $`i`$ and $`j[0,N]`$. To do so, the last three terms in Eq. (11) are written in a slightly different form,

$`W(\theta )^2^2\psi +2W(\theta )\stackrel{}{}W(\theta )\stackrel{}{}\psi `$ (33)
$`(_y\psi )_x\left[W(\theta )W^{}(\theta )\right]+(_x\psi )_y\left[W(\theta )W^{}(\theta )\right],`$ (34)
so that each new term can be computed with the help of the intermediate arrays. For instance the last term is discretized as

$$(_x\psi _{i,j})\left[(WW^{})_{i,j+1}(WW^{})_{i,j1}\right]/(2h).$$
(35)
## IV Equilibrium shapes

### A Comparison between analytical and phase-field results

We first check our code by computing the equilibrium shape of the crystal. In thermal equilibrium, $`u=\mathrm{\Delta }`$ everywhere, and the evolution equation for the phase field reduces to

$`\tau (\theta )_t\psi `$ $`=`$ $`\left[\psi +\lambda \mathrm{\Delta }(1\psi ^2)\right](1\psi ^2)+W(\theta )^2^2\psi `$ (38)
$`+2W(\theta )\stackrel{}{}W(\theta )\stackrel{}{}\psi (_y\psi )_x\left[W(\theta )W^{}(\theta )\right]`$
$`+(_x\psi )_y\left[W(\theta )W^{}(\theta )\right].`$
Initially, the solid is given a circular (or a square) shape and the undercooling is set to $`\mathrm{\Delta }=\mathrm{\Delta }_0`$. Time integration of Eq. (38) is performed and the interface velocity along the $`x`$ axis, $`V_x`$ is computed at regular time intervals. If $`V_x>0`$, $`\mathrm{\Delta }`$ is decreased and it is increased if $`V_x<0`$. The increment $`\delta \mathrm{\Delta }`$ is divided by a constant $`a>1`$ each time $`V_x`$ changes sign. This procedure is repeated as long as $`\delta \mathrm{\Delta }`$ is larger than some prescribed value, $`\delta \mathrm{\Delta }_{min}`$. This scheme is known to converge to an equilibrium state .

The analytical equilibrium shape is given by

$$\{\begin{array}{cc}\stackrel{~}{x}=x(\theta )/R_0=f_s(\theta )\mathrm{cos}\theta f_s^{}(\theta )\mathrm{sin}\theta ,\hfill & \\ \stackrel{~}{y}=y(\theta )/R_0=f_s(\theta )\mathrm{sin}\theta +f_s^{}(\theta )\mathrm{cos}\theta ,\hfill & \end{array}$$
(39)
where $`R_0=d_0/\mathrm{\Delta }`$ . Let us consider one eighth of this interface corresponding to $`\theta [\pi /4,\pi /2]`$ (Fig. 2). In the cusp region, $`\theta [\pi /2\theta _0,\pi /2]`$, one has

$$\{\begin{array}{cc}\stackrel{~}{x}(\theta )=B\mathrm{cos}\theta ,\hfill & \\ \stackrel{~}{y}(\theta )=(B\mathrm{sin}\theta A),\hfill & \end{array}$$
(40)
so that the interface is the circle of center $`(0,A)`$ and radius $`B`$. Since the smoothed cusps are very narrow in practice ($`\theta _01`$) this circle arc is very flat; it tends to the straight horizontal facet as $`\theta _00`$. The right end of this interface portion lies at

$$\{\begin{array}{cc}\stackrel{~}{x}_r=\stackrel{~}{x}(\pi /2\theta _0)=B\mathrm{sin}\theta _0,\hfill & \\ \stackrel{~}{y}_r=\stackrel{~}{y}(\pi /2\theta _0)=B\mathrm{cos}\theta _0A,\hfill & \end{array}$$
(41)
where the local slope is

$$\stackrel{~}{y}_r^{}(\stackrel{~}{x}_r)=\mathrm{tan}\theta _0.$$
(42)
On the other hand, for $`\theta [\pi /4,\pi /2\theta _0]`$ the crystal is bounded by the circle centered at point $`P_c(\delta ,\delta )`$ (i.e. $`x=y=\delta `$) of rescaled radius unity, for which

$$\stackrel{~}{R}=(\stackrel{~}{x}_t\stackrel{~}{x}_r)\sqrt{2}=1\sqrt{2}\mathrm{sin}\theta _0$$
(43)
is a good first-order approximation when $`\theta _01`$ (Fig. 2).

Note that, as $`\theta _00`$, one recovers the equilibrium shape for the sharp cusp, with a horizontal facet of length

$$\stackrel{~}{\mathrm{\Lambda }}=2\stackrel{~}{x}_r2\delta $$
(44)
matching tangentially a circle of rescaled radius unity.

As shown in Fig. 3, the whole crystal shape is well reproduced by the phase-field equilibrium code. A close examination of the numerical data points reveals that the imposed anisotropy is underestimated by about 0.2 percent in the numerics, independently of $`\theta _0`$.

### B Numerical estimates for the facet ends

The knowledge of the analytical equilibrium shape guides us to define a numerical procedure to extract estimates of the facet length $`\mathrm{\Lambda }`$ and corner radius $`R`$ from discrete interfaces obtained with the phase-field code. On the square mesh, the interface consists in a list of points $`P_i(x_i,y_i)`$. According to Eq. (42), the two facet ends are located at the points $`P_r`$ and $`P_l`$ where the absolute value of the local slope rapidly increases beyond $`\mathrm{tan}\theta _0`$ (Fig. 4). The method used to compute the two end points is thus to compare the derivative of the equilibrium curve $`_xy`$ with $`\mathrm{tan}\theta _0`$. Using centered differences,

$$_xy_i=\frac{y_{i+1}y_{i1}}{x_{i+1}x_{i1}},$$
(45)
gives an error of about 5 percent on both measures. On the other hand, using the one-sided approximations,

$$_xy_i|^{}=\frac{y_iy_{i1}}{x_ix_{i1}},$$
(46)
to compute $`P_r`$ and

$$_xy_i|^+=\frac{y_{i+1}y_i}{x_{i+1}x_i},$$
(47)
to compute $`P_l`$, estimates of $`\mathrm{\Lambda }`$ and $`R`$ are only 0.5 percent off the exact values (see Table I ). This last estimation procedure is thus more precise and it is readily extended to the case of growth shapes in what follows.

### C Sharp cusp limit

One may still wonder if the sharp cusp limit ($`\theta _00`$) is reacheable within our numerical approach. To answer this question, it is better to develop Eq. (38) in terms of the second partial derivatives of $`\psi `$. One gets then

$`\tau (\theta )_t\psi `$ $`=`$ $`\left[\psi +\lambda \mathrm{\Delta }(1\psi ^2)\right](1\psi ^2)`$ (49)
$`+W_0^2\left[C_{xx}_{xx}\psi +C_{xy}_{xy}\psi +C_{yy}_{yy}\psi \right],`$
with

$$C_{xx}=\{\begin{array}{cc}1+\delta (3|\mathrm{cos}\theta ||\mathrm{cos}^3\theta |+|\mathrm{sin}^3\theta |)+\delta ^2,\hfill & \\ A^2+B^2AB|\mathrm{cos}\theta |(2+\mathrm{sin}^2\theta ),\hfill & \\ B^2AB|\mathrm{sin}^3\theta |,\hfill & \end{array}$$
(50)
$$C_{xy}=\{\begin{array}{cc}2\delta (\sigma _s\mathrm{cos}^3\theta +\sigma _c\mathrm{sin}^3\theta )+2\delta ^2\sigma _s\sigma _c,\hfill & \\ 2AB\sigma _c\mathrm{sin}^3\theta ,\hfill & \\ 2AB\sigma _s\mathrm{cos}^3\theta ,\hfill & \end{array}$$
(51)
and

$$C_{yy}=\{\begin{array}{cc}1+\delta (3|\mathrm{sin}\theta ||\mathrm{sin}^3\theta |+|\mathrm{cos}^3\theta |)+\delta ^2,\hfill & \\ B^2AB|\mathrm{cos}^3\theta |,\hfill & \\ A^2+B^2AB|\mathrm{sin}\theta |(2+\mathrm{cos}^2\theta ),\hfill & \end{array}$$
(52)
for $`\theta _0<\theta <\frac{\pi }{2}\theta _0`$, $`\theta \theta _0`$, and $`\theta \frac{\pi }{2}\theta _0`$, respectively. The symbols $`\sigma _s`$ and $`\sigma _c`$ represent the signs ($`\pm 1`$) of $`\mathrm{sin}\theta `$ and $`\mathrm{cos}\theta `$.

The time step $`\mathrm{\Delta }t`$ must be sufficiently small to ensure convergence of the finite difference scheme, the most stringent constraint on $`\mathrm{\Delta }t`$ arising for $`\theta =\theta _0`$:

$$\mathrm{\Delta }t\frac{\tau _0}{W_0^2}\frac{h^2}{2}\left[\frac{f_s^2(\theta _0)}{C_{xx}(\theta _0)+C_{yy}(\theta _0)}\right].$$
(53)
As $`\theta _01`$, this approximates to

$$\mathrm{\Delta }t\frac{\tau _0}{W_0^2}\frac{h^2}{2}\left(\frac{1+\delta }{\delta }\right)\theta _0.$$
(54)
It is thus always possible to reduce $`\theta _0`$ provided that $`\mathrm{\Delta }t`$ is decreased in proportion. Although $`\mathrm{\Delta }t`$ depends on the cusp amplitude $`\delta `$, the time step can still be kept constant while maintaining numerical stability by imposing

$$\theta _0=K\frac{\delta }{1+\delta },$$
(55)
with $`K`$ some constant. In the numerical simulations, we tipically take $`K=\pi /100`$, so that $`\mathrm{\Delta }t/\tau _00.0025\mathrm{}`$. In practice, we rather use the discretization scheme described in section III than the fully developped version given in Eqs. (49-52). The former proves to be more precise and more stable than the latter, so that a larger time step, $`\mathrm{\Delta }t/\tau _0=0.008`$, can be used.

## V Faceted needle growth

We now turn to non-equilibrium growth shapes. The time evolution of a needle dendrite is illustated in Fig. 5 for a cusp amplitude $`\delta =1`$, diffusion coefficient $`D\tau _0/W_0^2=4`$ and undercooling $`\mathrm{\Delta }=0.55`$. The size of the simulation box is $`600W_0\times 600W_0`$. For reasons of symmetry, it is sufficient to grow one half of the dendrite (here $`yx`$). We checked that identical patterns are obtained when the whole quadrant is used. The initial conditions for the phase field are $`\psi =1`$ if both $`0x20W_0`$ and $`0y20W_0`$ and $`\psi =1`$ otherwise. This square shape is imposed to the solid germ in order to get a dendrite with a single tip along the $`x=y`$ direction. Conversely, rounded initial shapes, like a circle, always result into tip splitting along the growth direction, leading to more complicated patterns which will be discussed later. The temperature field $`u`$ is initially set to zero inside the solid germ and to $`\mathrm{\Delta }`$ outside. Reflecting conditions are imposed at all the domain boundaries

We first focus on the influence of the anisotropy coefficient, $`\delta `$, on the dynamics of the faceted needle dendrite. This coefficient is varied from 0.25 to 1.60 in the simulations. After a transient, the dendrite reaches a stationary state which is independent of the size of the initial germ. The tip then moves with a constant velocity $`V`$ and a rather large portion behind the tip adopts a stable shape. A good test of the thin-interface limit is to decrease independently the rescaled interface thickness $`W_0/d_0`$ (or equivalently the diffusion coefficient $`D\tau _0/W_0^2`$) and the grid spacing $`h/W_0`$ . Convergence is achieved when the scaled tip velocity $`\stackrel{~}{V}=Vd_0/D`$ is independent of both pameters. Table II shows that the results are reasonably well converged for the values $`D\tau _0/W_0^2=4`$ and $`h/W_0=0.4`$ used in our simulations.

### A Steady-state shape

The steady-state interface shape can be divided into three distinct parts: a nearly circular tip of radius $`R`$, a smooth “facet” of length $`\mathrm{\Lambda }`$, and a trailing rough tail (Fig. 6). Estimates of $`R=(x_tx_r)\sqrt{2}`$ and $`\mathrm{\Lambda }=x_rx_l`$ are obtained by using the numerical interpolation method described in section IV-B. To smooth out temporal fluctuations, time averages of both quantities are performed over the whole stationary regime. Obviously, the exact shape of the dendrite tip should depend on the smoothing angle $`\theta _0`$. When $`\theta _0`$ is sufficiently small this dependence is found to be linear, which allows us to unambiguously extrapolate $`R(\theta _0)`$ and $`\mathrm{\Lambda }(\theta _0)`$ to $`\theta _0=0`$ (Fig. 7). Up to an overall scale factor, the dendrite tip shape is very similar to the equilibrium shape, i.e. a quarter of a circle matching tangentially to two side facets.

The variations of $`\mathrm{\Lambda }(0)`$ and $`R(0)`$ with cusp amplitude are dispayed in Fig. 8. The results are well fitted by the two simple relations,

$$\mathrm{\Lambda }=\mathrm{\Lambda }_0\left[1+\frac{\delta }{\delta _\mathrm{\Lambda }}\right],$$
(56)
and

$$R=R_0\left[1+\frac{\delta _R}{\delta }\right].$$
(57)
The asymptotic shape of the trailing rough part away from the facet is well-fitted by a parabola that satisfies the Ivantsov relation

$$\mathrm{\Delta }=\sqrt{\pi }\sqrt{p}e^p\mathrm{erfc}(\sqrt{p})$$
(58)
with the Péclet number

$$p\frac{\rho V}{2D}.$$
(59)
Note that the reason why the asymptotic tail is a parabola whose tip radius is predicted exactly by the Ivantsov relation is the same as for non-faceted growth . Namely, in a frame that is stationary with respect to the melt, the trailing rough part grows as a planar interface whose position $`x(t)`$ perpendicular to the growth axis $`t^{1/2}`$ and the Péclet number $`px(t)(dx(t)/dt)/2D`$ obeys the same relation as Eq. (58), as shown independently of Ivantsov by Zener .

We checked in a few simulations that $`\rho V`$ is well predicted by the Ivantsov relation, as illustrated in Fig. 9 where we compare the phase-field interface to the Ivantsov parabola with the same tip velocity, i.e., the parabola with $`\rho =2Dp/V`$, where $`p`$ is the Ivantsov Péclet number and $`V`$ is the steady-state tip velocity in the phase-field simulation. Given this agreement, we use directly the formula $`\rho =2Dp/V`$ to calculate $`\rho `$ in all results reported thereafter. It should be emphasized the tip radius $`\rho `$ of this imaginary parabola that matches exactly the asymptotic needle crystal shape and the true tip radius $`R`$ are quite different as shown in the plot of $`R/\rho `$ versus $`\delta `$ in Fig. 10. In particular, $`R/\rho `$ is seen to decrease sharply with increasing $`\delta `$, while $`\mathrm{\Lambda }/\rho `$ approaches a constant. These plots reflect the fact that side facets become elongated and extend closer to the tip that becomes more pointed as $`\delta `$ increases.

### B Steady-state operating state

As in previous sharp-interface calculations , we define the dimentionless tip selection parameter

$$C\frac{4\rho ^2V}{Dd_0},$$
(60)
where $`V`$ is the tip velocity and $`\rho `$ is the tip radius of the parabola that exactly matches the asymptotic tail of the faceted needle crystal. Note that $`C=8/\sigma `$ where $`\sigma =2Dd_0/\rho ^2V`$ is the selection parameter that is also often used in dendrite growth studies.

Eliminating $`\rho `$ in favor of the Péclet number using the relation $`\rho =2Dp/V`$, we can rewrite $`C`$ in the form

$$C=\frac{16D}{Vd_0}p^2.$$
(61)
A plot of $`C`$ versus $`\delta `$ with $`V`$ corresponding to the steady-state tip velocity in the phase-field simulations and $`p`$ computed from the Ivantsov relation for the value of $`\mathrm{\Delta }`$ corresponding to the simulations is shown in Fig. 11. Note that $`C(\delta )`$ has a minimum that corresponds to a maximum of tip velocity as a function of $`\delta `$.

In experiments, the growth rate of a single material is usually studied as a function of undercooling Additional simulations are thus performed for other values of $`\mathrm{\Delta }`$. All the other parameters are kept constant, except for the total time which is increased to $`7000\tau _0`$ for $`\mathrm{\Delta }=0.45`$, $`3000\tau _0`$ for $`\mathrm{\Delta }=0.50`$, and reduced to $`1500\tau _0`$ for $`\mathrm{\Delta }=0.60`$. We do not extrapolate $`V`$ to $`\theta _00`$ here because the computing time becomes prohibitively large for $`\mathrm{\Delta }<0.55`$. However, from the previous results for $`\mathrm{\Delta }=0.55`$, we do not anticipate deviations in the tip velocity by more than a few percent (see Fig. 7). The plot of tip velocity versus facet length shown in Fig. 12 shows that our numerical results are consistent with the scaling law $`\mathrm{\Lambda }V^{0.5}`$ which was found experimentally for NH<sub>4</sub>Br needle crystals .

### C Transients

When growth is started from a circular germ, a double needle is systematically observed (Fig. 13), independently of the germ shape, for $`\delta 0.20`$. Initially, the two tips move away from each other. Our simulations show that these needles ultimately grow practically parallel to each other at very long times, when they are separated from each other by several diffusion lengths (Fig. 13).

At long times, the tip of the double needle becomes comparable in shape to that of the single needle dendrite (Fig. 14). Of course, the former is thinner on one side because the latent heat accumulates in the channel between the two tips. This similarity in shape suggests that the double needle pattern is governed by the same operating state as the needle dendrite. Verifying this point with a sufficient accuracy is not easy because of a still longer transient for the double needle dendrite. For this reason, we compare the patterns obtained for a larger undercooling, $`\mathrm{\Delta }=0.60`$, for which the simulations are still quantitative but the relaxation time much shorter than with $`\mathrm{\Delta }=0.55`$. After a time $`t=2000\tau _0`$, the tip velocity of the double needle pattern has not yet fully converged but it tends slowly to a limit near $`0.397W_0/\tau _0`$ (Fig. 15). This is precisely the value found for the single tip dendrite, which confirms that the two asymptotic states are identical. Thus, the double needle dendrite is merely just a tip-splitting evolution of the single tip dendrite. This splitting is often preferred by the system because it allows for a more efficient occupation of the available room.

Let us finally remark that the double needle found here is not a faceted version of the generic non-faceted growth structure that has been termed doublon . Doublons also are symmetry-broken growth shapes with a double tip. However, the two tips remain much closer than one diffusion length and are thus highly inter-dependent, which is not true for the double faceted needle dendrite observed in the present simulations.

## VI Comparison with sharp-interface analytical theory

A precise quantitative comparison of the phase-field results of the last section with previous sharp-interface calculations  cannot be carried out because these calculations considered the one-sided model with zero diffusivity in the solid-phase while we simulated here a symmetric model with equal diffusivities in solid and liquid. In addition, our simulations are for finite Péclet number while the sharp-interface calculations of Refs.  are restricted to the small Péclet number limit ($`p0`$). Nonetheless, the phase-field results reproduce several key features of faceted needle growth that were predicted in the sharp-interface studies:

1. The ratio $`\mathrm{\Lambda }/\rho `$ of the facet length to the tip radius of the parabola that matches the asymptotic rough tail increases with the cusp amplitude $`\delta `$ and saturates to a value of order unity for large $`\delta `$.
2. The tip selection constant $`C`$ has a minimum for a value of $`\delta `$ of order unity. The value of $`C`$ that corresponds to this minimum is large ($`100`$).
3. For large $`\delta `$, the needle shape consists of two long facets joined by a nearly circular tip that match on to trailing rough parts that become parabolic far from the tip. This is precisely the shape proposed by Adda Bedia and Hakim .

In the remainder of this section, we develop a simple analytical theory of faceted needle growth that is based on the approximate shape proposed by Adda Bedia and Hakim  for large $`\delta `$, where facets are joined by a small circular tip of radius $`R`$. We take this approximation further by assuming that (i) it remains valid for $`\delta `$ of order unity, and (ii) the trailing rough parts remain parabolic all the way to the points at which they match tangentially the facets, as opposed to being parabolic only asymptotically far from the tip . These simplifications allow us to obtain a simple physical picture of faceted needle growth as well as explicit analytical predictions of the shape parameters $`R`$, $`\mathrm{\Lambda }`$ and $`\rho `$ and the velocity $`V`$ for arbitrary Péclet number, without necessitating the numerical solution of an integral equation as in Refs. . As we shall see below, these predictions agree reasonably well with our phase-field results despite our simple parametrization of the steady-state shape.

The complete determination of the needle shape ($`R`$, $`\mathrm{\Lambda }`$ and $`\rho `$) and the growth velocity $`V`$ for fixed $`\mathrm{\Delta }`$ and $`\delta `$ requires four independent relations. The first is the Ivantsov relation (58) that fixes the product $`\rho V`$. A second relation between $`R`$, $`\mathrm{\Lambda }`$ and $`\rho `$ is simply obtained by imposing that the circular tip and parabolic tails match tangentially the front and trailing ends of the facets, respectively, which yields at the relation

$$\rho =(\mathrm{\Lambda }+R)/\sqrt{2}.$$
(62)
We can check this relation by comparing the tip velocity in the phase-field simulations with the velocity $`V=2Dp/\rho `$, where $`p`$ is the Péclet number predicted by the Ivantsov relation and $`\rho `$ is computed using Eq. (62) together with values of $`\mathrm{\Lambda }`$ and $`R`$ in the simulations. Fig. 16 shows that these two velocities are in reasonably good quantitative agreement. The systematic deviation is most likely due to the fact that the crystal tail is not exactly parabolic very close to the facet.

The two other relations needed to complete our theory are obtained by imposing the Gibbs-Thomson condition (3) at the tip of the needle crystal and by using the integral form of this condition on the facet (Eq. 5). For this, we need an expression for the undercooling along the interface that is obtained straightforwardly using the known boundary integral relation

$`u(x)=\mathrm{\Delta }+p{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dx^{}}{\pi }}\mathrm{exp}\left(p[z(x)z(x^{})]\right)`$ (63)
$`\times K_0\left(p\sqrt{(xx^{})^2+(z(x)z(x^{}))^2}\right),`$ (64)
where $`K_0`$ is the zeroth order modified Bessel function, $`p`$ is given by the Ivantsov relation (58), and $`z(x)`$ is the interface shape with length in unit of $`\rho `$. For a circular tip and parabolic tails matching tangentially side facets that make a $`45^o`$ angle with respect to the growth axis, the $`x`$ coordinates of the front and trailing ends of the facets are $`x=\pm R/\sqrt{2}\rho `$ and $`x=\pm 1`$, and

$`z=\left({\displaystyle \frac{R^2}{\rho ^2}}x^2\right)^{1/2}{\displaystyle \frac{\sqrt{2}R}{\rho }}+{\displaystyle \frac{1}{2}},0|x|{\displaystyle \frac{R}{\sqrt{2}\rho }},`$ (65)
$`z=|x|+{\displaystyle \frac{1}{2}},{\displaystyle \frac{R}{\sqrt{2}\rho }}|x|1,`$ (66)
$`z={\displaystyle \frac{x^2}{2}},|x|1.`$ (67)
These equations describe rather accurately the phase-field needle dendrites (Fig. 17). Since the isothermal Ivantsov parabola is an exact solution of the steady-state growth problem without capillarity, the integral relation (64) is exactly satisfied for $`z=x^2/2`$ and $`u=0`$. In contrast, as illustrated in Fig. 9, the tip of the faceted needle crystal protrudes ahead of the tip of this parabola and hence is undercooled by a finite amount $`u(0)`$. The Gibbs-Thomson condition (3) imposes a relation between this tip undercooling and the tip radius that is simply

$$u(0)=d_0/R.$$
(68)
In turn, Eq. (5) provides us with a relation between the average undercooling on the facet and the cusp amplitude $`\delta `$ that takes the form here

$$\sqrt{2}_{\frac{R}{\sqrt{2}\rho }}^1𝑑xu(x)=2\frac{d_0}{\rho }\delta ,$$
(69)
where $`x`$ is in units of $`\rho `$ as above. The four relations defined by Eqs. (58, 62, 68) and (69), together with the expression for $`u`$ on the interface defined by Eqs. (64-67), completely determine $`R`$, $`\mathrm{\Lambda }`$, $`\rho `$ and $`V`$. Since the product $`\rho V`$ is exactly predicted by the Ivantsov relation, we only need to compute three independent dimensionless combinations of the four above quantities to compare the predictions of the above theory with the phase-field results. The most meaningful dimensionless combinations are the selection parameter $`C`$ and the two ratios $`\mathrm{\Lambda }/\rho `$ and $`R/\rho `$. To compute those as a function of $`\delta `$, it is convenient to vary $`R/\rho `$ and compute $`C`$, $`\delta `$, and $`\mathrm{\Lambda }/\rho `$ using the following relations that are simple to deduce:

$`{\displaystyle \frac{1}{C}}={\displaystyle \frac{\sqrt{2}}{16p\delta }}{\displaystyle _{\frac{R}{\sqrt{2}\rho }}^1}𝑑xu(x),`$ (70)
$`\delta ={\displaystyle \frac{\rho }{\sqrt{2}u(0)R}}{\displaystyle _{\frac{R}{\sqrt{2}\rho }}^1}𝑑xu(x),`$ (71)
$`{\displaystyle \frac{\mathrm{\Lambda }}{\rho }}=\sqrt{2}\left(1{\displaystyle \frac{R}{\sqrt{2}\rho }}\right).`$ (72)
The above predictions are compared with the phase-field results in Figs. 10 and 11. The quantitative agreement between the two is reasonably good but not exact because the steady-state shape that we have assumed here differs slightly from the true shape (Fig. 17) that is only exactly parabolic asymptotically far from the facet.

To conclude, we note that the expression for $`u(x)`$ can be simplified in the small Péclet number limit where Eq. (64) reduces to

$$u(x)=p_{\mathrm{}}^+\mathrm{}\frac{dx^{}}{2\pi }\mathrm{ln}\left[\frac{(xx^{})^2+(z(x)z(x^{}))^2}{(xx^{})^2+(x^2/2x^2/2)^2}\right].$$
(73)
Therefore, both $`u(0)`$ and the integral of $`u(x)`$ along the facet are simply proportional to the Péclet number in this limit and hence $`C`$, $`R/\rho `$, and $`\mathrm{\Lambda }/\rho `$ are only functions of $`\delta `$. Furthermore, it is straightforward to deduce that $`\mathrm{\Lambda }/\rho \sqrt{2}`$, $`R/\rho 0`$, and $`C\delta `$ in the large $`\delta `$ limit. Hence this theory predicts the scaling law $`\mathrm{\Lambda }V^{1/2}`$ that is a simple consequence of the fact that $`\mathrm{\Lambda }/\rho `$ and $`C\rho ^2V`$ are constants in the small Péclet number limit.

## VII Concluding Remarks

In conclusion, we have shown that the phase-field approach can be successfully extended to model the solidification of faceted materials. Our approach, which consists of rounding the cusps in the $`\gamma `$-plot, converges well in the limit of sharp cusps for both equilibrium and non-equilibrium growth shapes with facets. Even though we have considered a simple form of the $`\gamma `$-plot, the method should be applicable to more complex $`\gamma `$-plots where the interface stiffness varies on the rough parts.

In addition, we have developed an approximate analytical theory of faceted needle growth that includes capillarity and assumes circular and parabolic forms for the front and trailing rough parts of steady-state needle crystals, respectively. This theory yields explicit predictions of tip velocity and facet length that are in good overall quantitative agreement with the phase-field results; this agreement is largely due to the fact that the needle shape assumed in the theory is a very good approximation to the actual growth shape observed in the phase-field simulations. Furthermore, this theory predicts that the scaling law $`\mathrm{\Lambda }V^{1/2}`$ observed both experimentally  and in the present phase-field simulations should hold in the small Péclet number limit, consistent with the results of previous theoretical studies .

In the present study, we have established a quantitative comparison between phase-field and sharp-interface results under the assumption that the variation of interface undercooling is dominated by capillarity. Using scaling arguments, Ben Amar and Pomeau  have concluded that, under the assumption that the growth of facets is dominated by a Franck-Read screw discolation mechanism, the kinetic undercooling of the interface should be negligible over an intermediate range of velocity where the $`\mathrm{\Lambda }V^{1/2}`$ scaling law has been observed . While growth conditions where kinetic effects are negligible may exist, we suspect that such effects will generally be important in the presence of facets, as widely believed. Therefore, a quantitative incorporation of facet kinetics in a phase-field model remains an important task for the future. The phase-field model as formulated thus far reproduces a linear relationship between the planar interface velocity and interface undercooling appropriate for a rough interface. On facets, however, standard growth mechanisms such as screw dislocations or ledge nucleation lead to a non-linear relationship between velocity and undercooling. These relationships could potentially be incorporated in the phase-field model by letting the kinetic relaxation time $`\tau `$ depend on temperature or supersaturation, in addition to orientation.

Another interesting future prospect is to model the directional solidification of alloys with faceted interfaces  by combining the present methodology to handle cusped $`\gamma `$-plots with a recent thin-interface phase-field formulation of alloy solidification .

## ACKNOWLEDGMENTS

A.K. wishes to thank the hospitality of the Laboratoire Matériaux et Microélectronique de Provence, Université d’Aix-marseille III where this work was initiated as well as the support of U.S. DOE Grant No. DE-FG02-92ER45471. We also thank Vincent Hakim and Klaus Kassner for valuable discussions. This research benefited from computer time allocation at the Advanced Scientific Computation Center at Northeastern University.

Tables: