Densest binary sphere packings and phase diagram : revisited

As discussed by Pickard and coworkers [33], randomly generated structures may contain spheres that are very close together, or spheres may be crowded in a direction of a short lattice vector instead of a long lattice vector. To avoid generating such an abnormal initial structure, we develop the piling up method to generate appropriate initial structures randomly. The method is based on an idea of stacking spheres randomly one by one on top of a randomly generated first layer, since any periodic structure can be understood as a multi-layered structure if it is extended in a direction perpendicular to a chosen base plane.

The piling up method consists of two steps. At first, an initial structure is generated randomly. Next, a multi-layered structure is constructed by expanding initial structures.

The generated multi-layered structure contains a large gap, which makes the packing fraction reduced. The large gap can be filled by repeating collision and repulsion between spheres under pressure. The steepest descent method with the hard-sphere potential defined as Eq. (10) can cause a repetition of collision and repulsion. If the constant repulsive force is large enough, the pressure and the repulsive force cannot be balanced. The repulsion is dominant when overlapping, leading to a structure without any overlap as a result of the optimization. Once the optimization reaches a structure without any overlap, the pressure is only the driving force to change the structure because of the characteristic feature of the two-body interaction defined by Eq. (10), and the optimization leads to a structure with overlap again. This means that the optimization cannot finish, but enable collision and repulsion between spheres to be repeated as long as we continue. The pressure is necessary to minimize the volume of the unit cell. In each steepest descent step, the extended coordinates $\bm{u}$ is changed to minimize the enthalpy per unit cell $H$ as

where $k_{2}$ is the steepest descent prefactor and the enthalpy is defined as $H\equiv E+PV$; $P$ is the pressure and $V$ is the volume of the unit cell. The operation is not aimed at converging the enthalpy to a local minimum but transforming significantly an initial structure to a dense structure with a large number of the repetition of collision and repulsion. In fact, we have confirmed that the linear minimization method does not work. A lot of iteration is essential for sufficient structural transformation. In the global optimization step, the steepest descent method is repeated several thousand times.

In the optimization, no structure stops transforming, because the forces and the pressure are never balanced. Any structure can transform into a dense structure. The effectiveness comes from the hard-sphere potential defined as Eq. (10).

As with the expansion process discussed before, the steepest descent prefactor $k_{2}$ is not set directly, but $\Delta\bm{u}$ is scaled so that the maximum displacement in the position vectors and the lattice vectors will be equal to the predetermined value. The purpose of the global optimization is only to fill the wasteful gaps, so it is the most important that how much the position vectors and lattice vectors are allowed to displace. Setting the maximum displacement has the advantage of allowing each sphere to move enough even when the number of spheres in the unit cell is large. Therefore, it makes sense to set $\Delta\bm{u}$ to an appropriate size so that the wasted gaps can be filled. During the global optimization, the maximum displacement is kept constant. The pressure $P$ is also kept constant as $P=0.1$. The other values are presented in the Sec. III.1.

Even after a large number of optimization steps, the overlap never converges to zero, because $z_{ij}^{\left(\bm{T}\right)}$ between neighboring spheres oscillate around zero. The overlap converges to zero with the local optimization, which we discuss the next.

To optimize a structure to a periodic packing of nonoverlapping spheres, all of the overlaps have to be converged to zero. In other words, all of the $z_{ij}^{\left(\bm{T}\right)}$ have to be converged to more than or equal to zero. To achieve the condition, we develop the iterative balance method. The iterative balance method optimizes the volume of a unit cell while keeping the overlap of hard spheres minimized. The method is developed based on the idea of repeating collision and repulsion among spheres under pressure while the maximum displacement in position vectors and lattice vectors is gradually decreased. The pressure $P$ is kept constant as $P=0.1$ through the iterative balance method.

In the iterative balance method, collision and repulsion between spheres are also repeated with the steepest descent method. As in the global optimization, the steepest descent prefactor $k_{2}$ is not set directly, but $\Delta\bm{u}$ is scaled so that the largest displacement in the position vectors and lattice vectors will be equal to a certain value. To converge all of $z_{ij}^{\left(\bm{T}\right)}$ to more than or equal to zero, the maximum displacement is gradually decreased for each step.

To illustrate the convergence of $z_{ij}^{\left(\bm{T}\right)}$, we consider the structural optimization in the one-dimensional case. Here, as shown in Fig. 1, we assume that two spheres interact with the discontinuous force defined to be

with

where $r$ is the distance between the two spheres. If the model is optimized with the steepest descent method with a constant prefactor, $z$ oscillates around $z=0$, and never reaches to $z=0$. On the other hand, if the displacement is gradually reduced, $z$ eventually converges to zero. The convergence behavior is shown in Fig. 1 together with the one-dimensional model and the discontinuous force defined by Eq. (15) as in the inset.

The idea works even for the three-dimensional case we are interested in. If the repulsive force is much larger than the applied pressure $P$, some of $z_{ij}^{\left(\bm{T}\right)}$ between neighboring spheres oscillate around zero when the steepest descent method with a constant prefactor is applied. Once we control the prefactor so that the maximum displacement can be gradually reduced, the structure eventually reaches a packing structure with almost zero overlaps. In this sense, the $z_{ij}^{\left(\bm{T}\right)}=0$ can be regarded as a balanced point. This is the reason why we call it the iterative balance method. The size of the maximum displacement is the same order as the maximum overlap which is defined as the maximum value in $-z_{ij}^{\left(\bm{T}\right)}$, as discussed in Sec. III.2. It should be clearly noted that the discontinuous sign change in forces produced by Eq. (10) is crucial to realize the balanced point. Other potentials such as the squared potential of $z$ do not provide the interesting feature.

The iterative balance method can find optimal distortion. Due to the pressure, the volume of a unit cell is decreased as much as possible. Therefore, the method makes as many $z_{ij}^{\left(\bm{T}\right)}$ converge to zero as possible, because the void between spheres causes an increase in the volume. Some of $z_{ij}^{\left(\bm{T}\right)}$ between neighboring spheres converge to zero. The convergence causes distortion in the structure. Accordingly, the generated packing structure has a high packing fraction, because the volume of the unit cell is minimized. The feature indicates that the iterative balance method can find the local maximum of packing fractions despite the fact that a large number of optimization steps is necessary to find the maximum packing fraction. The optimization parameters such as the optimization step number are presented in the Sec. III.1.

Just like the Torquato-Jiao sphere-packing algorithm [20], in principle the iterative balance method allows us to evaluate the packing fraction of the dense sphere packings while avoiding overlaps of spheres. The implementation is straightforward and the computational time for the structural optimization is very short as discussed later on, since it requires only simple calculations of forces and stress at each steepest decent step. We will demonstrate the efficiency and ability to find putative densest structures in Sec. III.

The (7-3) structure [22] is the tower-like DBSP with seven small and three large spheres in the unit cell. The formation process of the structure is shown in Fig. 2. Firstly, the tower-like structure is generated by stacking many spheres on top of the small bottom plane as shown in Fig. 2(a), and then the structure is expanded as shown in Fig. 2(b) with the steepest descent method. The bottom plane corresponds to the first layer generated by the piling up method. The final structure, shown in Fig. 2(c), also has a small bottom plane and is long in the vertical direction.

The (16-4) structure is the DBSP with sixteen small and four large spheres in the unit cell. The structure can be understood as a two-layered structure except for four small spheres. The formation process of the structure is shown in Fig. 3. For visibility of the figures, the bottom plane, which is the first layer generated by the piling up method, is arranged at the front side. Firstly, all the spheres are placed on the first layer as shown in Fig. 3(a), and then the structure is expanded as shown in Fig. 3(b). An overlap in the direction perpendicular to the first layer causes a distribution of spheres’ positions in the direction. The final structure, shown in Fig. 3(c), can also be regarded as a layer-by-layer structure with a small height.

If a packing structure sufficiently converges to a dense packing, few spheres’ positions displace significantly even after many optimization steps. Besides, repulsive forces occur only when two spheres overlap each other. Therefore, the computational cost can be reduced significantly by predetermining the neighboring spheres. Of course, we must reset the neighboring list periodically, because accidentally position vectors may displace largely. In our code, the neighboring list is reset once every 100 times in the global optimization step and once every 200 times in the local optimization step.

One may consider that it is futile to optimize all similar structures which converge to an isomorphic structure. Hereafter, the isomorphic is defined that two structures are equal when the distortion is corrected. The distortion is necessary to achieve a high packing fraction. In some cases, the number of local minima is huge due to a large number of distortion patterns and accordingly many structures are trapped at local minima. Therefore, optimizing all of the similar structures is necessary to determine the highest packing fraction. Hence, all of the generated structures are optimized with our method.

The efficiency of structural optimization depends on the choice of parameters used in the optimization, and even the reachable maximum packing fraction is varied depending on the parameters in case that many local minima are having competitive packing fractions. These parameters we have are listed below:

• •

  $\Delta_{\mathrm{expand,max}}$ : The maximum displacement in the expansion step

• •

  $\Delta_{\mathrm{global,max}}$ : The maximum displacement in the global optimization step

• •

  $N_{\mathrm{global}}$ : The maximum step number in the global optimization step

• •

  $\Delta_{\mathrm{local,max}}$ : The maximum displacement in the local optimization step

• •

  $N_{\mathrm{local}}$ : The maximum step number in the local optimization step

• •

  $d$ : The decreasing factor of the maximum displacement in the local optimization step

The default values used in our study are given below.

Here, we note that $\Delta_{\mathrm{global,max}}$ ($\Delta_{\mathrm{local,max}}$) is chosen randomly in the range of $0.03\leq\Delta_{\mathrm{global,max}}\leq 0.15$ ($0.02\leq\Delta_{\mathrm{local,max}}\leq 0.05$) for each structural optimization. The random choice of them enables us to obtain various local minima. The default parameters lead to a convergence that almost all the minimum value of $z_{ij}^{\left(\bm{T}\right)}$ becomes the same order as $-10^{-7}$ to $-10^{-6}$.

In the iterative balance method, the maximum displacement of the position vectors and the lattice vectors are gradually reduced. Figure 4 shows the maximum overlaps at each optimization step. $\Delta_{\mathrm{local,max}}$ is set to $0.0546375$. The maximum overlap is defined as the maximum value of $-z_{ij}^{\left(\bm{T}\right)}$. It is confirmed from Fig. 4 that the maximum overlap decreases exponentially as the maximum displacement is decreased exponentially, and those values are the same order. The relation can be understood from the fact that some of $z_{ij}^{\left(\bm{T}\right)}$ between neighboring spheres oscillate around zero.

In the radius ratio of $\alpha\leq\sqrt{2}-1$, the $\mathrm{XY_{n}}$ structures appear on the phase diagram. The $\mathrm{XY_{n}}$ structures are defined as packing structures in which large spheres X constitute the fcc densest structure and the small spheres penetrate into the tetrahedral and octahedral sites constituted by X. This definition is the same as that of Hopkins et al. [22]. The packing fractions of the $\mathrm{XY_{n}}$ structures can be calculated analytically as

In the radius ratio of $\alpha\leq\sqrt{2}-1$, as discussed by Hopkins and coworkers, six $\mathrm{XY_{n}}$ structures of $\mathrm{XY}$, $\mathrm{XY_{2}}$, $\mathrm{XY_{4}}$, $\mathrm{XY_{8}}$, $\mathrm{XY_{10}}$, and $\mathrm{XY_{11}}$ [22] appear on the phase diagram. An octahedral site of $\mathrm{XY}$ structure is occupied by one small sphere, an octahedral site of $\mathrm{XY_{2}}$ structure is occupied by two small spheres, an octahedral site of $\mathrm{XY_{4}}$ structure is occupied by a tetrahedron consisting of four small spheres, an octahedral site of $\mathrm{XY_{8}}$ is occupied by an cubic consisting of eight small spheres, a tetrahedral site of $\mathrm{XY_{10}}$ structure is occupied by one small sphere and an octahedral site of $\mathrm{XY_{10}}$ is occupied by a cubic consisting of eight small spheres, and a tetrahedral site of $\mathrm{XY_{11}}$ structure is occupied by one small sphere and an octahedral site of $\mathrm{XY_{11}}$ is occupied by a bcc structure consisting of nine small spheres, respectively.

It was relatively easy to calculate the packing fractions of the $\mathrm{XY_{n}}$ structures with our method because there is no distortion in the fcc structure constituted by large spheres. However, when the optimization parameters are set to the default values, in some cases an error of the packing fraction becomes the same order as $10^{-6}$ due to overlap which is also the same order as $10^{-6}$. As discussed in the next section, by optimizing again with changing the optimization parameters, the maximum overlap can be decreased to the same order as $10^{-13}$. In that case, the packing fractions agree by more than 10 decimal points with the analytical solution of Eq. (17).

As discussed in Sec. II.3, the iterative balance method can minimize the volume of a unit cell. Therefore, high-precision re-optimization enables us to find the highest packing fractions with high accuracy because of a small overlap and a minimized volume of a unit cell.

As discussed in Sec. III.3, the local optimization with default values presented in Sec. III.1 cannot reduce an overlap enough. The overlap causes an error in packing fraction after the six decimal points. The error can be reduced by a high-precision re-optimization. With an increase in the number of local optimization steps, the overlap can be decreased because the maximum displacement is decreased as the number of optimization steps is increased. As discussed in Sec. III.2, the maximum overlap and the size of the maximum displacement is the same order.

When the number of local minima is getting larger, more structures are trapped at local minima, so the global minimum might not be found during an exhaustive search. Therefore, it is necessary to re-optimize the putative densest structure many times in order to identify the global minimum. The optimization results depend on the initial structure and the optimization parameters such as $\Delta_{\mathrm{global,max}}$ and $\Delta_{\mathrm{local,max}}$. Therefore, slight fluctuations are given to the structure before re-optimization and the range of $\Delta_{\mathrm{global,max}}$ and $\Delta_{\mathrm{local,max}}$ is set to be large. When the number of local minima is small, almost all the fluctuated structures converge to the densest packing. On the other hand, when the number of local minima is large, we obtain many packing fractions corresponding to diverse distortion patterns. In some cases, tens of thousand re-optimizations are necessary to determine the maximum packing fraction. As discussed in the next subsection, the computational time for the structural optimization is very short, so the necessity of the several thousand re-optimization is not a serious problem.

The default parameters for re-optimization are given below.

If we use the default parameters, in almost all case, the minimum values of $z_{ij}^{\left(\bm{T}\right)}$ become the same order as $-10^{-13}$ to $-10^{-12}$.

The computational cost for calculating forces with the hard-sphere potential is very low because repulsive forces occur only when two spheres overlap each other. Besides, if a packing structure converges to a dense packing, the positions of spheres do not change significantly even after many optimization steps. Therefore, as discussed in Sec. II.5, the computational cost can be significantly reduced by pre-determining the neighboring spheres. In fact, if the number of spheres in the unit cell is about ten, it takes less than 0.1 seconds from generation to optimization using a single core of Intel(R) Xeon(R) CPU E5-1650 v4 @ 3.60GHz. Therefore, our method can generate a large number of packing structures. The efficiency enables us to find the densest packing from the vast coordination space.

In addition, since the number of optimization steps is set externally, the computational cost can be estimated as $O\left(N\right)$, where $N$ is the number of spheres. Hence, it is possible to conduct an exhaustive search for the long periodic densest packings.

In this subsection, we discuss the distribution and updating history of packing fractions during the exhaustive search, to show that our method can find the densest packings from the vast coordination space.

First, we conducted an exhaustive search for the mono-phase densest binary sphere packing at the radius ratio of $\alpha=0.445$ and the composition ratio of $x=14/19$. The unit cell contains 14 small spheres and 5 large spheres. In that case, the (14-5) structure with the packing fraction of $0.765259$ is the densest. The number of structures generated during the exhaustive search is 250,182. The distribution of the packing fractions is shown in Fig. 5(a). The number of structures having the packing fraction between 0.6875 and 0.6900 is 15959 and the largest. On the other hand, the number of structures having the packing fraction between 0.7650 and 0.7675 is only 58. The result indicates that our method can find the densest packings from the vast coordination space. Figure 6(a) shows the updating history of the highest packing fraction. The structure with a packing fraction of $0.765256$ is generated at the 9372nd steps, while the densest structure with the packing fraction of $0.765259$ is generated at the 50155th steps.

Similarly, we conducted an exhaustive search for the mono-phase densest binary sphere packing at the radius ratio of $\alpha=0.455$ and the composition ratio of $x=16/20$. The unit cell contains 16 small spheres and 4 large spheres. In that case, the (16-4) structure with the packing fraction of $0.759629$ is the densest. The number of structures generated during the exhaustive search is 351,124. The distribution of the packing fractions is shown in Fig. 5(b). The number of structures having the packing fraction between 0.6850 and 0.6875 is 24303 and the largest. On the other hand, the number of structures having the packing fraction between 0.7575 and 0.7600 is only 8. The result also indicates that our method can find the densest packings from the vast coordination space. Figure 6(b) shows the updating history of the packing fraction. The densest structure with the packing fraction of $0.765259$ is generated at the 51097th steps.

To determine the densest phase separation at each composition ratio $x$ and each radius ratio $\alpha$, firstly the mono-phase densest binary sphere packings have to be identified at each $\left(x,\alpha\right)$.

To identify the mono-phase densest binary sphere packings at each $\left(x,\alpha\right)$, up to one million structures are generated with SAMLAI. We terminate the exhaustive search when the highest packing fraction is not updated 200,000 or 300,000 times, and regard the structure with the highest packing fraction as putative densest packing.

The radius ratio $\alpha$ is changed by a step of 0.02 in the range of $0.20\leq\alpha\leq 0.64$ while the radius ratio is changed by a step of 0.005 in the radius range of $0.42\leq\alpha\leq 0.50$. In total, the 35 radius ratios are investigated in the construction of the phase diagram.

The number of spheres in the unit cell is set between 6 and 24 while the number of spheres in the unit cell is set between 12 and 32 in the radius range of $0.42\leq\alpha\leq 0.50$. Under the condition, all possible compositions are investigated with the constraint that the number of small spheres is equal to or larger than that of large spheres. In the former case, there are 138 compositions and in the latter case, there are 220 compositions, respectively. In general, the number of possible compositions is $n^{2}/4$ if the maximum number of spheres is set to $n$, where $n$ is an even number and the minimum number of spheres is set to 2.

The optimization parameters are set to the default as discussed in Sec. III.1.

Almost all of the generated structures have an overlap of about $10^{-7}$ to $10^{-6}$. Therefore, packing fractions may have an error on the sixth digit. In addition, the highest packing fraction may not be found due to a large number of local minima. Therefore, putative densest packings have to be re-optimized for determining the highest packing fractions with high accuracy. In some cases, tens of thousand steps for the re-optimization is necessary in order to find the highest packing fraction from the huge number of local minima. The re-optimization parameters are set to the defaults discussed in Sec. III.4 in most cases. The re-optimization updates the fourth digit and beyond of the packing fraction for cases that many local minima exist, while the packing fraction is unchanged by the re-optimization in most of the cases.

As Hopkins and coworkers have shown, the highest packing fractions for the binary system are achieved by phase separation into two or fewer densest packings [22]. We also proved it in a way that seems more intuitive to us. The proof is given in Appendix A. The phase separations with the highest packing fraction are determined from all possible phase separations.

We have discovered 12 putative densest packings and accordingly updated the phase diagram over the $x-\alpha$ plane, where $x$ is the relative concentration and $\alpha$ is the radius ratio of the small sphere relative to the large sphere. The phase diagram is shown in Fig. 7 for $0.20\leq\alpha\leq 0.40$, Fig. 8 for $0.420\leq\alpha\leq 0.500$, and Fig. 9 for $0.52\leq\alpha\leq 0.64$, respectively. On the phase diagrams, 24 DBSP are plotted. Some DBSP are not plotted, e.g., newly found (12-1) structure and the (9-4) structure, because they appear in a very narrow region. For the case of 12 or fewer spheres in the unit cell, our phase diagram is consistent with the HST phase diagram [22] with a small correction.

The $\mathrm{XY_{n}}$ structures are defined as DBSP in which the large spheres X constitute the fcc densest structure and the small spheres penetrate into the tetrahedral and octahedral sites constituted by X. If the $\mathrm{XY_{n}}$ structures are excluded, there are 21 putative DBSP. Most of those structures are named as ($m$-$n$) structure. A ($m$-$n$) structure contains $m$ small spheres and $n$ large spheres. Their packing fractions are shown in Tables 2, 3, and 4.

For binary systems, the highest packing fraction is generally achieved by phase separation into two or fewer densest packings. When the highest packing fraction is achieved by phase separation into a structure $\Phi$ and the fcc densest packing consisting of small spheres, only the symbol of structure $\Phi$ is plotted on the phase diagram. On the other hand, when the highest packing fraction is achieved by phase separation into two structures, those two symbols are plotted together on the phase diagram. For example, in the area of $0.30\leq\alpha\leq 0.34$ and $6/7<x$ the highest packing fraction is achieved by the phase separation into the (6-1) structure and the fcc densest structure consisting of small spheres, and in the area of $0.30\leq\alpha\leq 0.34$ and $1/2<x<6/7$ the highest packing fraction is achieved by the phase separation into the $\mathrm{XY}$ structure [22] and the (6-1) structure [22].

For $0.20\leq\alpha\leq 0.40$, the exhaustive search is conducted for 11 radius ratios of 0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, and 0.40. The number of spheres in the unit cell is set between 6 and 24. By the exhaustive search, we identify 6 new putative densest packings: the $\mathrm{XY_{12}}$ structure, the (22-1) structure, the (12-1) structure, the (20-1) structure, the (8-2) structure, and the (2-1) structure. The $\mathrm{XY_{12}}$ structure and the (22-1) structure appear on the phase diagram at $\alpha=0.20$ and $0.203$. The (12-1) structure appears on the phase diagram at $\alpha=0.203$; it can be understood as an extended $\mathrm{XY_{12}}$ structure. The (12-1) structure is not plotted on the phase diagram of Fig. 7, because it appears in a very narrow region. The (8-2) structure and the (20-1) structure appear on the phase diagram at $\alpha=0.26$; it is distorted $\mathrm{XY_{4}}$ structure. The (2-1) structure appears on the phase diagram at $\alpha=0.28$; it is distorted $\mathrm{XY_{2}}$ structure.

For $0.420\leq\alpha\leq 0.500$, the exhaustive search is conducted for 17 radius ratios of 0.420, 0.425, 0.430, 0.435, 0.440, 0.445, 0.450, 0.455, 0.460, 0.465, 0.470, 0.475, 0.480, 0.485, 0.490, 0.495, and 0.500. The number of spheres in the unit cell is set between 12 and 32. By the exhaustive search, we identify 5 new putative densest packings: the (14-5) structure, the (16-4) structure, the (8-4) structure, the (10-4) structure, and the (9-4) structure. The (14-5) structure appears on the phase diagram at the three radius ratios of $\alpha=0.440$, $0.445$, and $0.450$. The (16-4) structure appears on the phase diagram at the three radius ratios of $\alpha=0.450$, $0.455$, and $0.460$. The (8-4) structure appeares on the phase diagram at $\alpha=0.450$ and $0.465$ and it is isomorphic to the $\mathrm{HgBr_{2}}$ structure. The (10-4) structure appears on the phase diagram at $\alpha=0.480$ and $0.481$; it is isomorphic to the (5-2) structure. The (9-4) structure appears on the phase diagram at the three radius ratios of $\alpha=0.481$, $0.482$, and $0.483$. In addition, in contrast to the HST phase diagram, there is a phase separation into the (7-3) structure [22] and the $\mathrm{HgBr_{2}}$ structure at $\alpha=0.470$. Finally, the packing fractions of the $\mathrm{AuTe_{2}}$ structure [13, 22] are found to be consistent with those of the (4-2) structure [22]. Otherwise, the results are consistent with the HST phase diagram [22].

For $0.52\leq\alpha\leq 0.64$, the exhaustive search is conducted for 7 radius ratios of 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, and 0.64. The number of spheres in the unit cell is set between 6 and 24. By the exhaustive search, we identify one new putative densest packing: the (12-6) structure. It is isomorphic to the $\mathrm{AlB_{2}}$ structure [13, 22]. Otherwise, our phase diagram is consistent with the HST phase diagram [22].

In Secs. V.2, V.3, and V.4, the densest packings in each radius range are detailed.

For $0.20\leq\alpha\leq 0.40$, the exhaustive search is conducted for 11 radius ratios of 0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, and 0.40. In order to confirm that the (12-1) structure appears on the phase diagram, the $\mathrm{XY_{12}}$ structure, the $\mathrm{XY_{11}}$ structure, and the (22-1) structure are re-optimized at $\alpha=0.203$. The number of spheres in the unit cell is set between 6 and 24. As a result, we have identified the new six putative DBSP: the $\mathrm{XY_{12}}$ structure, the (22-1) structure, the (12-1) structure, the (20-1) structure, the (8-2) structure, and the (2-1) structure. Thereby, we update the phase diagram. For the case of 12 or fewer spheres in the unit cell, our phase diagram is consistent with that of the previous study [22] with a small correction. The putative DBSP already discovered are the $\mathrm{XY}$ structure [22], the $\mathrm{XY_{2}}$ structure [22], the $\mathrm{XY_{4}}$ structure [22], the $\mathrm{XY_{8}}$ structure [22], the $\mathrm{XY_{10}}$ structure [22], the $\mathrm{XY_{11}}$ structure [22], the (10-1) structure [22], the (11-1) structure [22], and the (6-1) structure [22]. On the phase diagram of Fig. 7, the $\mathrm{XY_{8}}$ structure, the $\mathrm{XY_{11}}$ structure, and the (11-1) are not plotted, because they appear in the narrow region.

In this subsection, firstly we present the small modifications of the HST phase diagram for the case of 12 or fewer spheres in the unit cell. Secondly, we present the newly found six putative DBSP. The four of them have extended unit cells compared to the previous study [22].

For $0.420\leq\alpha\leq 0.500$, the exhaustive search is conducted for 17 radius ratios of 0.420, 0.425, 0.430, 0.435, 0.440, 0.445, 0.450, 0.455, 0.460, 0.465, 0.470, 0.475, 0.480, 0.485, 0.490, 0.495, and 0.500. In order to confirm that the (9-4) structure appears on the phase diagram, the (2-2) structure, the (10-4) structure, the $\mathrm{AuTe_{2}}$ structure and the (9-4) structure are re-optimized at $\alpha=0.481$, $0.482$, and $0.483$. The number of spheres in the unit cell is set between 12 and 32. As a result, we have identified 5 new putative DBSP: the (14-5) structure, the (16-4) structure, the (8-4) structure, the (10-4) structure, and the (9-4) structure. Thereby, we update the phase diagram. For the case of 12 or fewer spheres in the unit cell, our phase diagram is consistent with that of the previous study [22] with a small correction. The putative DBSP already discovered are the (6-6) structure [22], the $\mathrm{HgBr_{2}}$ structure [13, 22], the (2-2) structure [22], the (5-2) structure [22], and the (7-3) structure [22] shown in Fig. 12, and the $\mathrm{AuTe_{2}}$ structure [13, 22] shown in Fig. 13.

In this subsection, firstly we present small modifications of the HST phase diagram for the case of 12 or fewer spheres in the unit cell. Secondly, we present the newly found five putative DBSP. The four of them have extended unit cells compared to the previous study [22].

For $0.52\leq\alpha\leq 0.64$, the exhaustive search is conducted for 7 radius ratios of 0.52, 0.54, 0.56, 0.58, 0.60, 0.62, and 0.64. The number of spheres in the unit cell is set between 6 and 24. As a result, we have identified new one putative DBSP, named the (12-6) structure. The structure have an extended unit cell compared to the previous study [22]. Thereby, we update the phase diagram. For the case of 12 or fewer spheres in the unit cell, our phase diagram is completely consistent with that of the previous study [22]. The putative DBSP already discovered are the $\mathrm{AuTe_{2}}$ structure [13, 22] shown in Fig. 13, and the $\mathrm{AlB_{2}}$ structure [13, 22] and $\mathrm{A_{3}B}$ structure [14, 22] shown in Fig. 15.

The (12-6) structure shown in Fig. 16 appears on the phase diagram at $\alpha=0.54$ and $0.56$. The packing fractions are $0.780466$ and $0.779098$, respectively. The structure can be understood as a distorted $\mathrm{AlB_{2}}$ structure as well as the $\mathrm{AuTe_{2}}$ structure. At these radius ratios, the $\mathrm{AuTe_{2}}$ structure and the $\mathrm{AlB_{2}}$ structure have the same packing fraction; the equivalence indicates that a large unit cell is necessary to realize an optimal distortion for higher packing fraction. At the three radius ratios of $\alpha=0.58$, $0.60$, and $0.62$, the (12-6) structure and the $\mathrm{AlB_{2}}$ structure have the same packing fraction; however, if the unit cell is expanded, a distorted $\mathrm{AlB_{2}}$ structure might appear on the phase diagram. The further investigation will be in future work.

The putative DBSP consist of local structures with high packing fraction. For example, $\mathrm{XY_{n}}$ structure is clearly made up of a combination of local structures with high packing fractions: the tetrahedron and octahedron constituted by large spheres. The same is true for the distorted $\mathrm{XY_{n}}$ structure such as the (2-1), (8-2), (10-1), (11-1), and (22-1) structures. The same is also true for the (6-6) structure in which small spheres penetrate into the octahedral sites in the distorted hcp densest structure constituted by large spheres. Similarly, the $\mathrm{AlB_{2}}$ structure, the $\mathrm{AuTe_{2}}$ structure, and the (12-6) structure consisting of the dense local structure in which a small sphere is embedded in the center of a triangular prism consisting of large spheres.

Some of the putative DBSP are complex structures, especially for $0.42\leq\alpha\leq 0.50$. Both the $\mathrm{HgBr_{2}}$ structure and the (7-3) structure consist of equilateral triangular prisms and parallelepiped hexahedrons. The (10-4) structure consists of parallelepipedal hexahedrons and cubics constituted by large spheres with five small spheres. The (16-4) structure consists of $\mathrm{AlB_{2}}$-type local structure and cubic frameworks constituted by large spheres with an octahedron made of small spheres. In the (14-5) structure, 14-oligomer structures constituted by small spheres are embedded in the gap among the large spheres. As discussed below, these local structures also appear in the mono-phase densest binary sphere packings that have a lower packing fractions than the densest phase separations. The appearance shows that those local structures are dense.

When space is filled with two kinds of spheres under periodic boundary conditions, the highest packing fraction is achieved by phase separation into two or fewer packing structures. However, in the below discussion, we prohibit a phase separation. Under the restriction, a lot of mono-phase densest sphere packings are obtained.

Mono-phase densest binary sphere packings also consist of some of the dense local structures. The (14-8) structure is shown in Fig. 17(a). The structure is the densest at $\left(\alpha,x\right)=\left(0.44,14/22\right)$. The unit cell contains 14 small spheres and 8 large spheres. The structure contains 14-oligomer structures, which is the same local structure in the (14-5) structure. The (14-8) structure contains a wasteful gaps compared to the (14-5) structure, but the presence of the (14-5)-type local structure suggests that the local structure is dense. The (11-6) structure is shown in Fig. 17(b). The structure is the densest at $\left(\alpha,x\right)=\left(0.48,11/17\right)$. The unit cell contains 11 small sphere and 6 large spheres. The structure contains the (10-4)-type local structure. The (11-6) structure contains a wasteful gap compared to the (10-4) structure, but the presence of the (10-4)-type local structure indicates that the local structure is dense.

Some of the mono-phase densest sphere packings are long-period structures; those structures are phase-separated into a few local structures in the unit cell. We discuss the two examples as such cases. Firstly, the (9-5) structure is shown in Fig. 18. The structure is the densest at $\left(\alpha,x\right)=\left(0.46,9/14\right)$. The unit cell contains 9 small spheres and 5 large spheres. It consists of the $\mathrm{HgBr_{2}}$-type triangular prism phase and the (6-6)-type stacking phase. The appearance of the $\mathrm{HgBr_{2}}$-type local structure indicates that the local structure is dense at $\alpha=0.46$. Secondly, the (16-6) structure is shown in Fig. 19. The structure is the densest at $\left(\alpha,x\right)=\left(0.56,16/22\right)$. The unit cell contains 16 small spheres and 6 large spheres. It consists of the $\mathrm{AlB_{2}}$-type triangular prism phase and the remainder of small spheres. The appearance of the $\mathrm{AlB_{2}}$-type local structure indicates that the local structure is dense at $\alpha=0.56$. The structural feature of the (16-6) structure indicates that the highest packing fraction at $\left(\alpha,x\right)=\left(0.56,16/22\right)$ is achieved by phase separation into the $\mathrm{AlB_{2}}$ structure and the fcc densest structure constituted by small spheres. In fact, our phase diagram shows that the highest packing fraction at $\left(\alpha,x\right)=\left(0.56,16/22\right)$ is achieved by the phase separation into the (12-6) structure and fcc densest packing of small spheres.

Those results show that all of the mono-phase densest binary sphere packings are made of a few local structures with high packing fractions. The concept of dense local structures are very similar to the compact packing; the concept is defined for disc packings [36]. The compact packing corresponds to a structure in which all the circles are in contact with the perimeter and the gaps consist only of curvilinear triangles. It has been proved analytically that there are 164 ternary compact packings [36]. A local structure with a high packing fraction in the three-dimensional DBSP may be regarded as an extended concept of the compact packing for three dimensions.

Sometimes the number of local minima becomes large. In this section, we discuss some representative examples, where the complexity of local minima can be understood.

Firstly, we analytically compute the packing fraction of the $\mathrm{HgBr_{2}}$ structure [13, 22] shown in Fig. 12(b). Hereafter, the large radius is set to 1 and the small radius is set to $r$. The structure contains equilateral triangular prisms constituted by large spheres; the length of one side of the equilateral triangle is 2 and the height of the sides of the triangular prism is $h$. Small spheres are placed in the center of each side of the triangular prism and one of them is inserted inside until it touches the other two small spheres. We name the side as $K$. In that case, $h$ can be calculated as

The $\mathrm{HgBr_{2}}$ structure also contains parallelepipedal hexahedrons. The bottom surface of the hexahedron is side $K$. One of the large spheres contacts the four large spheres which constitute the side $K$; here we assume that the large sphere does not overlap with the inserted small sphere. In that case, the height of parallelepiped hexahedron $l$ can be calculated as

The distance $d$ between the large sphere and small sphere can be calculated as

The $\mathrm{HgBr_{2}}$ structure consists of two equilateral triangular prism and one parallelepiped hexahedron, so the packing fraction can be calculated as

If we set the small sphere radius $r$ as $0.46$, $d$ and $\phi$ are calculated as

On the other hand, the packing fraction calculated by SAMLAI is 0.759011. This small difference comes from the slight gap between small spheres; $d$ is slightly larger than 1.46. The gap allows the two small spheres placed in the center of the sides of the equilateral triangular prism to penetrate slightly into the interior of the triangular prism; the distortion allows a slight increase in the packing fraction.

Secondly, the $\mathrm{HgBr_{2}}$ structure and the (7-3) structure [22] shown in Fig. 12(e) are very similar. These two structures consist of equilateral triangular prisms and parallelepiped hexahedrons. The packing fraction is increased by small spheres penetrating slightly into the interior of triangular prisms. The parallelepiped hexahedron consisting of large spheres is necessary for periodic penetrating. In the region where the radius is getting larger, the (7-3) structure appears on the phase diagram instead of the $\mathrm{HgBr_{2}}$ structure. This is because the (7-3) structure has twice as many triangular prisms as $\mathrm{HgBr_{2}}$ structure; the connection of many triangular prisms is necessary for small spheres to penetrate into gaps in triangular prisms without making wasted void between large spheres.

Thirdly, the (10-1) structure [22] appears on the phase diagram at the two radius ratios of $\alpha=0.240$ and $0.260$. Table 2 shows the packing fractions of the (10-1) structure at several radius ratios. The (10-1) structure is an expanded or a distorted $\mathrm{XY_{10}}$ structure so that small spheres can penetrate into the gap constituted by large spheres. The distortion pattern is varied according to the radius ratio. At most radius ratios, the (10-1) structure is distorted, but at $\alpha=0.239$, $0.240$, and $0.242$, it becomes an enlarged $\mathrm{XY_{10}}$ structure whose symmetry is kept $Fm\overline{3}m$. We are surprised to find that the highly symmetric structure appears again at those radius ratios. This suggests that a delicate balance among the numerous distortion patterns determines which one has the highest packing fraction.

As we have discussed, the packing fraction is increased by distortions. Determinating the maximum packing fraction becomes more difficult when the number of distortion pattern is larger. Besides, the unit cell is also related to the difficulty in determining the highest packing fraction.

Most of the densest packings are distorted; however, if the distortions are adjusted, some structures become highly symmetric and their unit cell could be reduced. For example, the unit cell of the $\mathrm{AuTe_{2}}$ structure can be reduced to that of the $\mathrm{AlB_{2}}$ structure [13, 22] if the distortion of the $\mathrm{AuTe_{2}}$ structure is adjusted. The unit cell of the $\mathrm{AlB_{2}}$ structure and the $\mathrm{AuTe_{2}}$ structure are shown in Fig. 20(a) and Fig. 20(b), respectively. In the $\mathrm{AuTe_{2}}$ structure, the two large spheres in the unit cell are not in the same plane, so the highest packing fraction of the $\mathrm{AlB_{2}}$ structure cannot be larger than that of the $\mathrm{AuTe_{2}}$ structure. Furthermore, even if the unit cell of the $\mathrm{AlB_{2}}$ structure is expanded, the highest packing fraction of the $\mathrm{AuTe_{2}}$ structure is not always achievable. This is a representative example that unit cell places a limit on the maximum packing fraction that can be reached. Hence, several isomorphic structures consisting of different unit cells need to be optimized to determine the maximum packing fraction. This is one of the main reasons why it is difficult to determine the maximum packing fraction.

As discussed in the previous section, it is difficult to determine the maximum packing fraction. However, not only almost all the packing fractions calculated by our method are the same as those shown in the previous study [22] by more than three decimal points, but also some of the packing fractions have successfully been updated. The results indicate that the iterative balance method can find the maximum packing fraction. Besides, the computational time for the structural optimization is very short as discussed in Sec. III.5. The effectiveness enables us to find the densest packing from the vast coordination space. Finally, as discussed in Sec. II.3, the implementation is straightforward since it requires only simple calculations of forces and stress at each steepest decent step. Those validity indicates that the iterative balance method might also be useful for the optimization of higher-dimensional sphere packings.

As discussed in Sec. III.6, our method can create diverse packing structures, so our method can discover a wide variety of mono-phase densest packings from the vast coordination space, as discussed in Sec. V.5. The validity indicates that the piling up method is an effective way to generate initial structures for predicting densest packings.

Three-dimensional data of the DBSP and the SAMLAI package, in which our methods are implemented, are available on the website [37]. The distribution of the program package and the source codes follow the practice of the GNU General Public License version 3 (GPLv3).

In some cases, atoms in crystals are approximated as spheres: In ion-bonded materials such as NaCl [4], atoms are often spherically symmetrical because atoms have closed-shell structures due to the charge transfer. In intermetallic compounds such as AgCu [5], atoms are also sometimes spherically symmetrical due to the bonds formed by electrons populating in s-orbitals. In materials under high pressure, distances between atoms become so close that the directional orientation of the bond is weakened due to the strong repulsive force by Pauli’s exclusion principle [6]. Therefore, we can assume that many crystals can be understood as densest sphere packings [4]. Then, we investigated the correspondence of DBSP with the crystal structures with reference to the space groups. The Spglib [38] is used for determining the space group of a densest sphere packing. The distortions of DBSP are corrected by the Spglib. As a result, we have succeeded in finding many crystals corresponding to DBSP, e.g., $\mathrm{LaH_{10}}$ [34] and $\mathrm{SrGe_{2-\delta}}$ [35] synthesized under high pressure, $\mathrm{YH_{6}}$ [32] predicted theoretically under high pressure, $\mathrm{NaCl}$; their crystal structures correspond to the (4-2) structure [22], the $\mathrm{XY_{10}}$ structure [22], the (6-1) structure [22], and the $\mathrm{XY}$ structure [22], respectively. Besides, the crystal structure of $\mathrm{UB_{4}}$ corresponds to the (16-4) structure. All of the correspondence we find are shown in Table 1. These correspondences indicate that the densest sphere packings can be used effectively as structural prototypes for searching complex crystal structures, especially for high-pressure phases.

On the other hand, we could not find the correspondence between crystals and the 15 densest binary sphere packings: $\mathrm{XY_{4}}$, $\mathrm{XY_{8}}$, $\mathrm{XY_{11}}$, $\mathrm{XY_{12}}$, (8-2), (11-1), (12-1), (20-1), (22-1), (14-5), (8-4), $\mathrm{HgBr_{2}}$, (7-3), (10-4), and (9-4). However, we consider that these structure can be realized by crystals, e.g., the (20-1) structure and the (22-1) structure may be realized by hydrides. As discussed in Sec. V.5, the DBSP consist of a few local structures with high packing fractions. This is a different principle compared to crystals, but it is apparently true that the volume of a crystal structure is decreased as much as possible under high pressure, so many of the densest packings may correspond to crystals especially for high-pressure phases.