The materials with the specific physical properties, in particular with negative refractive index play a significant role in the process of improving the radiation performances of the different IC and radio electronic devices. Such materials are used widely for improving the directivity characteristics of microstrip antennas of the different types, microwave filters and field transformers. There are different approaches to form the specific properties of the medium (material) by embedding into it a series of particles; this leads to forming the physical properties of resulting material that are different from those inherent to properties of the initial one. The theoretical prediction of existence of such materials was made in the pioneering work  and starting from this time such materials were designed by the various recipes. As early as the eighties of the last century, these materials got the name “chiral” and began to be used in various areas of antenna technology , the manufacture of electronic devices  , and telecommunications equipment .
The usual approach for solving the scattering problems for the acoustic and electromagnetic applications foresees the use of the integral equation method with the subsequent adaptation of different projection methods. In the main, such combination is used for the bodies and domains with the coordinate surfaces. The approaches noted above were used effectively for solving the problems of elastostatics and elastodynamics  , electrodynamics , as well as for a series of applied problems related to engineering and industry . A comprehensive description of these approaches was made in monograph , where the problems of both a theoretical nature and examples of specific applications in the various fields of science and technology are considered. But the disadvantage of the above methods is that they can be effectively used for the bodied/domains with the coordinate surfaces if the analytical approach is used or a huge computational resource is needed if the problem is solved for the complex geometry. Therefore, application of the asymptotic method that allows obtaining an explicit solution to the considered problem, which does not depend on the form of scattering object, is very promising.
The goal of this paper is to propose the numerical approach (based on  and ) for creating the material with the specific refractive index including its negative values. The approach foresees reducing an explicit asymptotic solution to the respective acoustic scattering problem, and the explicit formula for the resulting refractive index based on the asymptotic solution to acoustic wave scattering problem on a set of a big number of embedded particles of small size.
The paper is organized as follows. Section 2 is devoted to statement of diffraction problem and outline of application limits of geometrical and physical parameters of the material under consideration. The analytical form of solution will be derived in Section 3; and the numerical aspects of solving the respective system of linear algebraic equations (SLAE) will be presented. In Section 4, the explicit formula for the refractive index of resulting inhomogeneous material will be derived. The numerical results, related to exactness of asymptotic solution to the initial diffraction problem and properties of obtained refractive index will be presented in Section 5. A short conclusion finalizes the discussed topic under consideration.
2. Statement of Scattering Problem
A combination of both the asymptotic method and numerical simulation is used to solve the problem of creating a material with specific scattering characteristics, particularly, with a given refractive index. The initial diffraction problem is solved by assumption: , and , a is the characteristic size of the particle, d is the distance between adjacent particles, is the wave number.
An asymptotic solution to the problem of scattering on many particles by assumptions: , and was obtained in ; M is the total number of particles contained in a given domain .
The impedance boundary conditions
are prescribed on the boundary of m-th particle, where is boundary impedance, , ; is arbitrary function, continuous in ; , and , where , and .
The incident field satisfies Helmholtz equation in the whole , by this the scattered field satisfies the radiation conditions. We assume that a small particle is sphere of radius a with center in point , .
The full field satisfies equation
in the domain (2)
and boundary conditions
in , where (3)
and full field is
is solution to problem (2)-(4) at , (namely, when D contains no particles), is incident field, and field satisfies the radiation conditions.
Let belong to , and is arbitrary subdomain of D; is number particles in determined by
where function is prescribed and continuous in domain D.
It was substantiated in  that there exists some specific field (limiting field), which satisfies the next condition
and solution to the initial diffraction problem (2)-(4) can be sought for from the equation
where is the Green function for Helmholtz Equation (2) for the case of absence the particles. This fact allows us to use the approximate solution instead of exact solution and to obtain an explicit formula for refractive index of the constructed inhomogeneous material.
3. The Closed Form Solution to Scattering Problem
In order to derive the explicit formula for approximate field, we introduce the concept of limiting (effective) field . In paper  it was proved that the exact solution to problem (2)-(4) can be presented in form
Despite the fact that the last presentation contains an unknown function in the integrand, in contrast to formula (7), where all functions in the integrands are known, it is used to obtain an approximate solution to the original scattering problem. For this goal, we define the effective field , which acts on the m-th particle as:
and the next relation for the neighboring points is valid . We present the exact formula (8) in form
where the values are
Using the known relation for function from , and the asymptotic representation for values , we obtain the next formula for
at for (12)
The values are defined by the asymptotic formula
and the asymptotic formula for function is
Using last two formulas, we obtain the asymptotic representation of the effective field in the vicinity of particles
which is valid in the domains , where .
In order to calculate the values of the effective field everywhere using formula (15), we should know the values . They can be easy obtained as the solutions to the following system of linear algebraic equations (SLAE):
where , and . The matrix of SLAE (16) is diagonally dominant; therefore it is convenient for numerical solving. It was proven in  that this SLAE has unique solution for sufficiently small a.
In order to justify the exactness of solution to SLAE (16), which is used for determination of the effective field, we derive some different SLAE, corresponding to limiting Equation (7). Let us divide the domain D, where the small particles are located, into an union of the small non-intersecting cubes with centers in points ; the diameter of such cubes can be chosen as . Because the limited quantity of cubes cannot give whole D, we consider their smallest fragmentation that contains D, and define values in the cubes, which do not belong to domain D.
In order to determine the solution to Equation (7), we apply for it the collocation method proposed in  and applied successfully to solving the similar SLAR in  . In accordance with this method, we obtain such SLAE
where P is number of cubes that form a partition of D, is a center of p-th cube, is its volume. Since the value of d is small, then diameter can be of an order larger than distance d between particles. Since , then solving SLAE (17) is much easier than solving SLAE (16) in terms of the number of calculations.
As a result, we have two different SLAE (16) and (17). Solving both the SLAE, we can compare their solutions and to evaluate the area of accuracy of asymptotic solution (15). This evaluation has also the practical importance because it allows to find the optimal parameters of the domain D, which provide the possibility to create the refractive index that is the closest to the desired one.
4. Refractive Index of Resulting Material
The explicit formula (7) for the effective field opens the way to determination of the refractive index of the obtained material. It is important from the practical point of view how the calculated refractive index differs from that is obtained from the theoretical assumptions. We confine here by the real refractive index and formulate the constructive algorithm to obtain desired refractive index. It consists of three steps.
Step 1. Using known and unknown , we calculate function .
Step 2. Using the relation we determine
Equation (18) for two unknown functions and has infinite number of solutions , for which the conditions are fulfilled.
In this connection, the solution to (18) we determine as:
Calculation of and by (19) finalizes Step 2 of our procedure.
Step 3 is completely constructive and its goal is the following
• to create on the small particle of radius a the necessary impedance ;
• to embed the particles that satisfies the properties (19) into domain D with the initial properties.
Let the domain D consists of a set of small nonintersecting cubes with center in the points , and we place in each cube the quantity
the small spheres with radius a in center with point . Here value defines the nearest integer to . Let us distribute the balls with the distance of , and prescribe the boundary impedance of spheres with the value , where function is calculated by (19). It was proven in  that the modeled material, which is obtained by such embedding of small particles into domain D, has the desired refractive index , and its error tends to zero at . This theoretical result justifies the numerical procedure, which is applied to calculating the values of refractive index for the resulting material and determination of the relative error of new refractive index. Numerical calculations performed allow us to study in more details the role of the domain D parameters ( ) in the process of forming the resulting refractive index of new inhomogeneous material.
The application of the above algorithm was considered in  for the case of complex function , the above algorithm can be applied if material is lossless.
5. Numerical Modeling
5.1. Checking the Applicability of Asymptotic Solution
5.1.1. Exactness of Solution of the Limiting Equation (7)
The computational checking for determination of the exactness of solution to (7) foresees carrying out the calculations with a set of different problem’s parameters. We calculate the absolute and relative errors in the process of growth of number of the collocation points. The dependence of error of the parameter p ( ), where P is total amount of the small domains in D for , , and for the different values of function are shown in Figure 1 and Figure 2. The solution, which corresponds to ( ), is considered as benchmark one.
The relative error of solution to (7) is equal to 1.05% and 0.053% for the real and imaginary parts if (53 collocation points), it diminishes to 0.72% and
Figure 1. The relative error versus parameter p, .
Figure 2. The relative error versus parameter p, .
0.023% respectively at (63 collocation points), and to values 0.29% and 0.018% at (83 collocation points), (see Figure 1). This error is less than 0.009% for the real part of solution if; it tends to zero if value p grows. The error depends of the values of function too; it diminishes if the imaginary part of decreases (see Figure 2); the error of real part of solution if is equal to 0.009%, and error of imaginary part is thousandths of a percent. The obtained results confirm that calculation of the values of approximate field can be carried out with the high enough accuracy, and this accuracy is attained in a wide range of the geometrical and physical parameters of the material under investigation.
The results of computations show that the relative error depends of the parameter k to a large extent. In Figure 3 and Figure 4, the error is shown at and respectively,. One can see that the error is of one order larger at. The maximal error (if) at is less on 27% than those for.
5.1.2. Comparison of Solutions to Limiting Equation (7) and Asymptotic SLAE (16)
In the previous Subsection, we consider the solution to SLAE (17) with as benchmark solution to Equation (7). The maximal value of relative error for this p does not exceed 0.009% for wide range of the problem’s parameters.
The numerical calculations are presented for the different sizes of D and different functions. The obtained results for the small values of m are shown in Table 1 for, , and (linear sizes of parameters a, d, and in the text and tables below are measured in cm). The values of, which are received by formula (21), when the expected number of particles is changed by M. For this case, the radius of particle is determined as
The values of in the third row correspond to the optimal values of radius a, which guarantees the minimal error for module of solution to Equation (7) and system (16). In the fourth row, the values of distance d between particles are shown. The error’s maximal value is attained at; the error diminishes, if m grows. The difference between the optimal value and estimated value of particle’s radius depends on the parameter m to a large extent. So, this difference is equal to 26.09%, 5.15%, and 4.61% at, , and respectively. The distance d between particles decreases if m grows because
Figure 3. The relative error versus parameter p,.
Figure 4. The relative error versus parameter p,.
Table 1. The optimal parameters of domain D at small m,.
diameter of domain D remains constant. The relative error for solution to SLAE (16) is big enough but it does not exceed 1% already at.
The numerical results for large m with the same initial data are shown in Table 2. If to compare the difference between the parameters and, one can see that the difference between them is small enough, namely difference is in fourth character after the dot, although the relative difference is equal to 3.95% at for example. Similarly to the results presented in Table 1, the distance d between particles diminishes if m grows and it is an order of magnitude higher than the values of a. The minimal error of solution is attained at (the total number M of particles is equal to 27.4625 × 104), and it is equal to 0.20%.
Table 3 contains the comparable results for with the same set of initial data. One can see that the relative error diminishes if the number M of particles increases (one should note that the relative error depends of parameters a and too). This error tends to relative error of solution to Equation (7) if the value of m becomes larger than 80, namely. If to compare the results presented here with those are given in Table 1 and Table 2, one can
Table 2. The optimal parameters of domain D at large m,.
Table 3. The optimal parameters of domain D at large m,.
see that the relative difference between values and is small in absolute values, and the relative error for solution to SLAE (16) is less by an order of magnitude. This testifies that the value of function greatly affects the accuracy of the results.
5.1.3. Investigation of Difference between Solutions to SLAE (16) and (17)
The comparison of solutions to SLAE (16) and (17) was carried out for the different values of a at the different p and m. The relative error of SLAE (16) diminishes if p increases while m remains constant. For example, if p increases on 50%, then the relative error diminishes on 11.7% (if and,).
The difference of solutions to SLAE (16) and (17) in the real part, imaginary part, and module is shown in Figure 5 and Figure 6 if and. By this, the difference of real parts does not exceed 3.9% at, it is less than 3.3% at, and it is less than 1.85% if, , and. Respectively, this difference is less than 0.075% if, , and, (m is the same).
The numerical results, obtained for a wide range of d, show that its optimal value exists; and starting from this value, the deviation of solutions begins to increase again. Such optimal values of d are presented in Table 4 for the different constant. The results obtained testify at the optimal distance d between particles increases if the number of particles grows. For the small number of particles, the optimal distance is of the same order that a, for the set number of particles (, namely), this distance is of order larger. The optimal
Figure 5. Dependence of difference of the solution’s components of distance d between particles,.
Figure 6. Dependence of difference of the solution’s components of distance d between particles,.
Table 4. Optimal values of d for the different constant.
distance d between particles depends also of the value of radius a; if a increases, this distance grows, so at the difference between optimal values of d is equal to 21.13%. The similar relation remains and for other values of N.
The values of the minimal and maximal errors, which are attained for the optimal d, are shown in Table 5.
Table 5. The relative error of solution to SLAE (16) in % (min/max) for the optimal d.
The obtained results allow to conclude that the optimal value of d diminishes slower, when function grows, what is more this decreasing is more significant for the smaller a. For those, who are interested in a more detailed study of accuracy of the obtained asymptotic results related to solving the limiting Equation (7) and the corresponding SLAE, author refers to monograph .
5.2. Modeling the Material with the Desired Refractive Index
Numerical calculations are conducted for the case if. For the simplicity, we consider the case when a given domain D consists of the same subdomains. This limitation is not essential for numerical modeling. Numerical calculations were performed for the case when, and, D is some cube with side and particles are placed uniformly in domain D (the relative error of the solution to system (16) does not exceed 0.01% for this P). Let the initial domain D be the material with the initial refractive index, and the desired refractive index be, namely constant. Then the values can be calculated by formula (20). On the other hand, we can choose the number m such that is closest to number. It is easy to see that the corresponding for such M is calculated by such formula
that is, the obtained value of the refractive index differs on. To obtain the minimum difference, we choose two numbers and that satisfy the inequality where and. Therefore, if we have the value for the fixed a, we can obtain the numbers and, and to calculate also the nearest to the values by the formula (22).
The dependence of the maximal relative error for the calculated values of on radius a of a particle is shown in Figure 7 for for complex function (in Figures 7-9, the solid and dashed lines correspond to the imaginary part and real part).
The obtained results testify that the relative error considerably depends of the parameters, , and. This error is smallest if the value of numbers and are much close to value. The error has periodic character that is defined by the properties of functions and by the values of parameters and. The mean of error in the period grows if a increases. The comparable results are shown in Figure 8 and Figure 9 at and respectively.
Figure 7. The maximal relative error of the modeled refractive index,.
Figure 8. The maximal relative error of the modeled refractive index,.
Figure 9. The maximal relative error of the calculated refractive index,.
The minimal value of error 0.79% is attained for and it is equal to 0.51% if, and 0.26% if for.
The uniform placement of particles in domain D is simplest from the engineering point of view. Using the data, given in Figures 7-9, we can evaluate the number M of particles, which is necessary to obtain the refractive index more closest to the desired one (at given parameter of domain D). The respective results are shown for in Figure 10. The values are shown in the y axis. The solid, dashed and dot-dashed lines correspond to values respectively. The knowledge about the optimal number of particles in the domain D is the subsequent step to creating a material with the given refractive index.
The data shown in Figure 10 testify that the optimal number of particles diminishes if their radius grows. The estimation determines the distance d between particles. This distance differs from that is determined by uniform placement of particles in D. For example, for, it is equal to 0.1361, and the calculated d is equal to 0.119 and 0.159 for and respectively. The computations show that the relative difference between these two values of d is almost proportional to the relative error of the refractive index.
Since this value d does not depend on size of D in accordance with estimation, it can be applied as the additional parameter for optimization while choosing the alternative values of m. On the other hand, we can evaluate the number of particles in D by formula (20). Having, we can easily calculate the quantity M of particles if they are distributed uniformly in D. The distance d between particles is calculated easily too if is prescribed.
The analysis of formula (22) for calculation of the resulting refractive index shows that its value does not depend of the radius a of particles, but this
Figure 10. The optimal value of m versus radius a of particles for the different.
parameter has significant role when determination of error of the desired and calculated refractive index. The results, presented in Figures 7-9 demonstrate that there are optimal values of a, which provide the minimal value of relative error. For example, the minimal value of this error for (see Figure 7) is equal to 0.46%, 0.63%, and 0.86% for, , and respectively; at the same time the maximal value of error is equal to 12.16%, 14.12%, and 19.25% for, , and respectively. This means that size of particles has a decisive role at the creating the refractive index closest to the desired one. One can see from the results, presented in Figures 7-9 that the relative error of the created refractive index diminishes if the radius a of particles decreases (this is especially evident in the values of the maximum error) that confirms the theoretical prediction obtained in .
The opportunity to obtain the negative values of resulting refractive index clearly follows from (22), namely one can choose such parameters M, , and function that the first term will be larger than the second one, and the resulting refractive index will be negative.
An asymptotic approach has been developed to solve the problem of acoustic scattering on a set of small size particles (bodies) placed in a homogeneous material. The scattering problem is reduced to solving a corresponding SLAE whose dimension is equal to the number of particles. The solutions of this system are used in the formula of explicit representation of the scattered field. Numerical calculations are performed that determine the accuracy of the obtained solution depending on the physical parameters of the problem.
The obtained numerical results demonstrate the possibility of applying the proposed technique to create materials with specified acoustic properties, in particular the refractive index. A constructive algorithm for modeling the material with the desired refractive index is proposed.
The results of numerical modeling open up the possibility of engineering solutions for practical applications. For example, uniform placement of particles is the easiest way to engineering design, and the answer to how many particles should be placed in a given domain is given by the numerical simulation results.
The engineering problems regarding the placement of a large number of small particles into a given domain D and creating on their surface of the necessary impedance require the separate technological solutions.
Author thanks to Prof. Alexander Ramm from Kansas State University, USA, for the continuous interest for the topic under investigation what made it possible to obtain a series of new physical results.
 Ogier, R., Fang, Y.M. and Svedendahl, M. (2015) Near-Complete Photon Spins Electivity in a Metasurface of Anisotropic Plasmonic Antennas. Physical Review X, 5, 041019-1-041019-8.
 Yang, Y., Da Costa, R.C., Fuchter, M.J. and Campbell, A.J. (2013) Circularly Polarized Light Detection by a Chiral Organic Semiconductor Transistor. Nature Photonics, 7, 634-638.
 Lachat, J.C. and Watson, J.O. (1976) Effective Numerical Treatment of Boundary Integral Equations: A Formulation for Three-Dimensional Elastostatics. International Journal for Numerical Methods in Engineering, 10, 991-1005.
 Bouchon, M. and Sánchez-Sesma, F. (2007) Boundary Integral Equations and Boundary Elements Methods in Elastodynamics. Advances in Geophysics, 48, 157-189.
 Chapko, R. and Johansson, T. (2018) A Boundary Integral Equation Method for Numerical Solution of Parabolic and Hyperbolic Cauchy Problems. Applied Numerical Mathematics, 129, 104-119.
 Constanda, C., Doty, D. and Hamill, W. (2016) Boundary Integral Equation Methods and Numerical Solutions. Springer International Publishing, Bern, Switzerland.
 Ramm, A.G. (2013) Electromagnetic Wave Scattering by Small Impedance Particles of an Arbitrary Shape. Journal of Applied Mathematics and Computing, 43, 427-444.
 Ramm, A.G. and Andriychuk, M.I. (2012) Scattering of Electromagnetic Waves by Many Thin Cylinders: Theory and Computational Modeling. Optics Communications, 285, 4019-4026.
 Andriychuk, M.I. (2018) Solving the Problem of Electromagnetic Wave Scattering on Small Impedance Particle by Integral Equation Method. Progress in Electromagnetics Research C, 81, 211-223.
 Andriychuk, M.I. and Ramm, A.G. (2010) Scattering by Many Small Particles and Creating Materials with a Desired Refraction Coefficient. International Journal of Computing Science and Mathematics, 3, 102-121.
 Andriychuk, M.I. (2019) Antenna Synthesis through the Characteristics of Desired Amplitude. Cambridge Scholars Publishing, Newcastle, UK.