Optimization of the composition in a composite material for microelectronics application using the Ising model

Tailored material is necessary in many industrial applications since material properties directly determine the characteristics of components. However, the conventional trial and error approach is costly and time-consuming. Therefore, materials informatics is expected to overcome these drawbacks. Here, we show a new materials informatics approach applying the Ising model for solving discrete combinatorial optimization problems. In this study, the composition of the composite, aimed at developing a heat sink with three necessary properties: high thermal dissipation, attachability to Si, and a low weight, is optimized. We formulate an energy function equation concerning three objective terms with regard to the thermal conductivity, thermal expansion and specific gravity, with the composition variable and two constrained terms with a quadratic unconstrained binary optimization style equivalent to the Ising model and calculated by a simulated annealing algorithm. The composite properties of the composition selected from ten constituents are verified by the empirical mixture rule of the composite. As a result, an optimized composition with high thermal conductivity, thermal expansion close to that of Si, and a low specific gravity is acquired.

www.nature.com/scientificreports/ the optimization problem can be attained by finding the lowest-energy state (ground state) of the Hamiltonian (energy function) in the Ising model using the adiabatic quantum computing approach 10 . In this paper, we introduce a new materials informatics approach in which the optimum composition in a material can be clarified using the Ising model which can be adopted by a quantum annealing system as well as a simulated annealing system.

Results and discussion
Composition optimization using Ising machine. Composite materials we examine in this paper are widely used in many applications, from aerospace to biomedical applications and microelectronics [11][12][13] . The composition of a composite material depends on the application because the requirement of the material varies with the characteristics of the final component product. To date, many studies regarding the empirical mixture rule of composite materials have been performed in metal-based composites as well as ceramic-based and resinbased composite materials 14 . Therefore, once the composition is fixed, the properties of the composite can be predicted by calculating the mixture rule. Also, Fuzzy Preference Selection Method (f-PSI) method and so forth are well-known as the method for selecting the optimal composite material from the given several composite candidates 15 . However, in the case of inverse problem, the best composition from among many constituent candidates for obtaining the designated properties, with all combinations obtained by changing the composition and constituent, must be calculated until the designed properties are obtained. Thus, since the more the constituent candidates and composition scale are increased, the more the combination is drastically increased, much time is required to calculate all combinations using the classical computer with the von Neumann architecture. For example, in the case of n = 25 cities in the traveling salesman problem, of which the combination is 3.1 × 10 23 (n!/2n) 16 , it is reported that it takes more than 1.3 billion years to complete all calculations even on a supercomputer, while it is solved in less than 1 min on a machine using heuristics optimization: the Ising model 17 . Therefore, to find the optimum composition of a composite material, we apply the same step in this study. Our approach consists of the following five steps.
Step 1: Extraction of problems.
Step 2: Conversion of combinatorial optimization problem.
Step 3: Formulation of Ising model expression.
Step 4: Conversion of expression into QUBO.
Step 5: Computation of QUBO for the optimal solution using Ising machine.
The detailed content of every step is described later. We choose a heat sink and heat spreader for a high-speed computer as a target component for the materials informatics approach using the Ising model. The heat sink and heat spreader are attached to the backside of the Si chip, establishing a CPU for the high-speed computer [18][19][20] , as shown in Fig. 1. Under operation, more than 50 W/chip, depending on the computer, is generated on the Si surface 21 . Large amounts of heat should be dissipated from the Si chip efficiently using a heat sink and a heat spreader component to properly operate the circuit on the active side of the Si. Therefore, the following three requirements must be satisfied by the material for a heat sink and heat spreader 22-24 . 1. High thermal conductivity for proper heat dissipation 2. A thermal expansion coefficient close to that of Si for maintaining the attachability with Si 3. A low specific gravity for reducing the mechanical load on the fragile Si and flip-chip bonding structure www.nature.com/scientificreports/ To meet the above three requirements simultaneously, we formulate an energy function using the quadratic unconstrained binary optimization (QUBO) format equivalent to the Ising model. The general equation of the Ising model regarding the Hamiltonian H is as follows 25,26 : where σ i represents an input variable satisfied as σ ∊ {− 1, + 1}. J ij is a (two-body) interaction given parameter, and h i is referred to as a given magnetic field as a (one-body) parameter. In practice, the model is converted to the equivalent QUBO style with bits of q ∊ {0, 1} instead of σ ∊ {− 1, + 1}. QUBO is obtained only by converting a variable in the Ising model [σ = 2q − 1 or q = (σ + 1)/2].
To adopt the method for solving the discrete optimization problem regarding the composition of a composite material, we introduce variable bit xij ∊ {0, 1} with a material composition (j) of material constituent (i). The values for i (m = 10) and j (n = 100) are applied. Then, m × n = 1000 xij bits are introduced in this solution. All bits are allocated as shown in Table 1.
Ten materials (m = 10) shown in Table 2 are selected as the constituent of the composite in this study. Two constrained terms must be introduced to properly map the bit pattern shown in Table 1.
One constraint term F requires that there be at most one bit that is 1 in a certain material constituent column i so that the bit represents a composition of the constituent. This constraint term F is represented by the following formula:  10 100 x 1,100 x 2,100 x 3,100 x **,100 x 9,100 x 10,100 A 100 www.nature.com/scientificreports/ Another constraint term G forces the sum of the compositions of the selected components to be m = 100. This constraint term G is represented by the following formula: In this study, three objective functions are necessary since three properties are required for the composite material. The first requirement is higher thermal conductivity. The objective function term of thermal conductivity E TC is expressed as follows: where O i represents the thermal conductivity of constituent i.
The second requirement is a thermal expansion close to that of Si: 3.6 ppm/℃. The objective function term of the thermal expansion E TE is expressed as follows: where P i represents the thermal expansion of constituent i.
The third requirement is a lower specific gravity. The objective function term of the specific gravity E SG is expressed as follows: where Q i represents the specific gravity of constituent i.
The total energy is formulated by adding three objective terms and two constrained terms in the following way: where α, β, γ, δ and ε are hyperparameters that adjust the balance between the constraint terms and the objective functions.
The optimum composition, which satisfies all three requirements simultaneously, is calculated by determining the bit pattern matrix in Table 1 so that the total E is minimized with the simulated annealing process of the Ising model machine 27 . The calculation is performed by one of the Ising machines, i.e., the digital annealer (DA) 28 , for which application-specific CMOS hardware is designed to solve fully connected QUBO problems, using a simulated annealing (SA) algorithm, with a massively parallel architecture to optimize the energy of the Ising model by a Markov chain Monte Carlo (MCMC) search 29,30 . The architecture processes 1024 parallel bits that are fully connectable through 16-bit weights.
As a result of the simulated annealing process with the Ising machine, the following bits are chosen as the optimized answer: x 2, 38 , x 5, 3 , x 6, 2 , x 7, 9 , x 9, 19 , x 10,29 . One of the optimum compositions in this study is 38Cu-3BeO-2AlN-9h-BN-19c-BN-29Diamond (percent by volume), depending on the hyperparameters α, β, γ, δ and ε. This is the final answer to the inverse problem solved by the Ising model machine. The role of hyperparameter is determining the function weight of every term in total energy equation and normalizing the energy value of every term. In this case, we set α: 0.001, β: 0.01, γ: 0.01, δ: 80 and ε: 100, respectively. The total calculation time of simulated annealing process was less than 10 s, like the 25 cities traveling salesman problem described above 16 . Verification of composite material properties using mixture rule. We verify the properties of the optimized composition of a composite material using the mixture rule 31 . The mixture rules for the thermal conductivity κ of composites, normally applied to composite materials, are shown below. The first equation is the mixing rule where the heat flux direction and the constituent materials are arranged in parallel, and the second equation is where the heat flux direction and constituent materials are arranged perpendicularly. The isotropic composites of the type in which constituent particles are distributed in a matrix match the empirical logarithmic law well κ i is the thermal conductivity of constituent material i, V i is the volume fraction of constituent material i.
In the mixture rules for thermal expansion coefficient α of the isotropic composites, the following Turner equation is well known www.nature.com/scientificreports/ α i is the thermal expansion coefficient of constituent material i, V i is the volume fraction of constituent material i, K i is the Young's modulus of constituent material i. The specific gravity ρ of a composite material can be calculated by the following equation if the secondary phases are not formed between constituent materials in a composite: ρ i is the specific gravity of constituent material i, V i is the volume fraction of constituent material i. Table 3 shows the calculated properties of the optimized composite material obtained by the mixture rule. The composite material with its composition optimized by the Ising model satisfies all three demanded properties, as shown in Fig. 2.
The composite material, which includes diamond, will be able tobe produced by conventional 32,33 and modern 3D-printing material processing 34,35 .
As described above, the target of this study is the material for the heat sink and heat spreader of a high-speed computer driven by a silicon active chip. Recently, similar types of heat sinks and heat spreaders have been needed for power electronics applications. For example, the power amplifier of the base station of the mobile terminal for 5G and post 5G generates much higher heat than that of a high-speed computer 36 . In that application, by adopting the thermal expansion of a material of the active device, i.e., SiC or GaN, the same energy function of the Ising model can be utilized.
Furthermore, other additional requirements such as the cost of material, eco-friendliness, and so forth, can be covered by adding the designated objective function term, though only three objective functions regarding thermal conductivity, thermal expansion and specific gravity are specified in this study. Moreover, this composition optimization approach is used for the dielectric composite material in the same manner, since the dielectric constant can be verified using the mixture rule as follows. www.nature.com/scientificreports/ The first equation is the parallel mixing rule, where the composite constituent material is aligned parallel to the electric field, and the second equation is the series rule, where the composite constituent material is aligned in series with the electric field. The isotropic composite constituent material is aligned randomly and empirically obeys the logarithmic law 37 ε is the dielectric constant of the composite, ε i is the dielectric constant of constituent material i, V i is the volume fraction of constituent material i.

Conclusions
We demonstrated that a new materials informatics approach that incorporates the Ising model is valuable for the optimization of the composition in a material. At this time, we use the Ising machine with fully connected 1024 variables in size. Since the variables are enhanced to 8192 by improving pairing application-specific CMOS hardware recently, much larger-scale materials informatic problem can be solved. Also, this approach can be applied to the quantum annealing machine because the QUBO format we studied is compatible with that machine.
This new materials informatics approach will become a more powerful tool as the development of quantum annealers and modified Ising machine enable more massive and faster calculations.

Methods
In this paper, no materials were used since a theorical approach was discussed. The methodology used was described throughout the paper.

Data availability
Authors can confirm that all relevant data are included in the paper.