Nonlinearity is one of the fundamental phenomena in nature and human society [1, 2]. By means of nonlinear mechanisms applied in artificial complex systems, various devices or concepts such as transistors, lasers, and artificial intelligence were created, leading to the dramatic revolution in modern science and technology [3–6]. On the other hand, nonlinearity in macro-scale heat transport systems is lacking study in both theory and applications [7, 8], although its counterpart at micro or nano scale has been a significant topic of phononics in the last two decades [9–11]. Comparable to the coupling-induced inharmonic interaction in phonon transfer [12], nonlinearity in macro-scale heat diffusion is mainly reflected in the intrinsic response to external fields, for example, thermal conductivity or capacity is varying with temperature [13, 14]. The absence of macroscopic phenomenological theory (under the Fourier’s law) makes it difficult to handle nonlinear parts in heat conduction and limits the regulation or management of heat in industrial engineering and daily life.

Thermal metamaterials have flourished as a promising scheme for manipulating heat since the proposal of transformation thermotics [15–19]. The range and sensitivity of accessible thermal conductivities are thus extended to the level far beyond natural materials [20–23]. More recently, various design methods have been generalized from linear to nonlinear systems, in which thermal parameters are temperature dependent [24–27], and lead to a broad category of smart or self-adapting thermal metadevices [28, 29]. In particular, composites of artificial architectures can result in a larger value on the coefficient of an effective nonlinear term than building-unit materials, which is usually called nonlinearity enhancement [8]. This effect in random [30], periodic [31], and core-shell [32] structures have been proposed. However, besides the enhanced parameter ranges, anisotropy in functionality is also a crucial benefit of utilizing artificial architectures [33–35], allowing several regulating abilities integrated into a single installation. But the effective thermal conductivities considered in the above works [30–32] are all in scalar form, i.e., isotropic. When composites have relatively strong anisotropy in configuration, such as ellipses inclusions, effective thermal conductivities will be anisotropic. Then the nonlinearity enhancement effect may be directional, depending on the intensity of anisotropy in composite media.

In this work, we aim at designing directional nonlinearity enhancement in composite thermal media. Directionality can be induced by elliptical particles embedded in homogeneous media with identical orientations. We build a two-dimensional theoretical model for deducing an analytical relation between directional enhancement and corresponding influencing factors, including the inclusive area fraction, shape factor, and ratio of linear part in intrinsic thermal conductivities, and demonstrate total-factor analyses with numerical methods. Finite-element simulations verify the designed model and give a visualized range and level of directional nonlinearity enhancement. Considering the universality of ellipses for mimicking a number of geometrical configurations such as clavae or circles by tuning shape factors, the proposed basic model may be extended to other anisotropic systems and inspire a broad category of multifunctional or Janus nonlinear thermal metadevices.

Let us consider a two-dimensional composite model in which a large number of ellipse inclusions are randomly distributed in a host matrix with identical orientation, see Figure 1. We use subscripts i, m, and e to indicate the parameters of inclusion, matrix, and effective medium, respectively. Then κ _i (T) and κ _m (T) are denoted to the intrinsic temperature-dependent thermal conductivities of two constituents, and κ _e (T) is the effective thermal conductivity of the composite. For simplification without loss of generality, intrinsic thermal conductivities are set to be composed of two parts, namely the linear and nonlinear components. They can be written in as κ i T = κ i 0 + χ i T α (1) and κ m T = κ m 0 + χ m T β . (2) The first terms on the right hand of the above two equations are linear parts (constants), while the second terms are nonlinear parts. χ is the nonlinear coefficient and T represents the local temperature. α and β can be assigned as arbitrary real numbers. It is noted that we consider the weak nonlinearity effect here, which is common for most solid crystals within the room temperature range. So κ _i0 ≫ χ _i T ^α and κ _m0 ≫ χ _m T ^β should be satisfied. Naturally, the Taylor expansion technique is adapted for deriving the analytical form of κ _e (T) from κ _i (T) and κ _m (T). Executing Taylor expansion on κ _e (T) by regarding χ _i T ^α and χ _m T ^β as small quantities, it can be expected to retain κ e T = κ e 0 + η 1 χ i T α + η 2 χ m T β + O χ i T α + O χ m T β , (3) where η ₁ χ _i T ^α and η ₂ χ _m T ^β are first-order expansion terms, and O(χ _i T ^α ) and O(χ _m T ^β ) are their higher-order expansion terms. By comparing first-order nonlinear terms of the effective medium and its components, we can obtain nonlinear gain coefficients. In our following analyses, we focus on the first-order nonlinear term, higher-order nonlinear terms are ignored. Eq. 3 implies that η ₁ and η ₂ are dimensionless nonlinear gain coefficients resulted from composite effects. It can be employed to evaluate the level of nonlinearity enhancement.

To proceed, we refer to the analytical form describing the effective thermal conductivity of the composite in the linear case, which is also known as the generalized Maxwell-Garnett (M&G) equation [36]. Under the case that the long axes of ellipse inclusions are along x axis, the nonlinear thermal conductivity in the x direction is expressed as κ e x T = κ m T + f κ i T − κ m T κ m T g x κ i T − κ m T + κ m T 1 − f g x κ i T − κ m T g x κ i T − κ m T + κ m T , (4) where f is the area ratio of inclusions to the whole media and g _x is the major-axis shape factor of ellipse inclusions, which is defined exactly in Ref. [36]. Similarly, the thermal conductivity in the y direction can be obtained by replacing x with y in Eq. 4. We can see the shape factor g is the key for constructing anisotropic effective thermal conductivities because g _x and g _y are different for κ _ex (T) and κ _ey (T), respectively. So we can naturally consider that g will induce divergent nonlinearity enhancement effects in different directions. Next, we will deduce the detailed form of effective nonlinear modulation coefficients. For defining them explicitly, we consider two simplified cases. One is that nonlinear inclusions embedded in linear matrix, the other is that linear inclusions embedded in nonlinear matrix.

To check our theory, we perform finite-element simulations solving the heat conduction equation with the commercial software COMSOL Multiphysics (www.comsol.com). The whole composite media is constructed as a square with a side length of 10 cm. A total of 400 ellipse particles are randomly embedded in the matrix. Each particle has a semi-axis length a in the x direction (see Figure 1), and the thermal bias ΔT = 1 K is also applied in this direction with the hot (cold) source at 301 K (300 K). For different cases, we take κ _i0 = 4 W m⁻¹ K⁻¹ and change the ratio κ _i0/κ _m0. To generate a weak nonlinearity, the nonlinear coefficient for the nonlinear material is set as 10^–4 W m⁻¹ K⁻². The effective nonlinear coefficient is calculated by comparing the effective conductivities when the thermal bias exists or is absent. In addition, the temperature in the effective nonlinear thermal conductivity is set as the average temperature over the inclusions, which is approximately equal to 300.5 K.

First, to show the (directional) nonlinearity enhancement by nonlinear inclusions and linear matrix, we give the simulated η with the theoretical results in Figure 4. According to Figure 2, we take g = 0.8 here to see the effect of nonlinearity enhancement in Figure 4A. For plots in different colors, we take κ _i0/κ _m0 = 0.1, 1/3, 1, 3, and 10, respectively. The data for the scatter plot are the average value of simulations using three different random position sets for the inclusions. The trend of the scatter plots basically agrees with the theoretical value. When κ _i0/κ _m0 = 0.1, we can see obvious nonlinearity enhancement when f ≥ 0.2. In addition, η is greater than the theoretical value when f is not small. This deviation comes from the overly simplistic assumption that the M&G theory only considers the dipole effect. In Figure 4B, we compare η with different g values while κ _i0/κ _m0 is fixed. The plots of g = 0.9 and g = 0.1 (or g = 0.7 and g = 0.3) tell that the directional nonlinearity enhancement do exists. Then, we give similar results for linear inclusions and nonlinear matrix in Figure 5. We can see the condition of κ _i0/κ _m0 for nonlinearity enhancement is different from the case in Figure 4. Anyway, we can observe (directional) nonlinearity enhancement as well.

Hereto, we have demonstrated the directional nonlinearity enhancement in anisotropic thermal media. Different from the previously-studied isotropic structures, non-circular configurations introduce directionality into effective nonlinear thermal conductivity, inducing different effective nonlinear coefficients in orthometric directions. In particular, the case that nonlinear inclusions embedded in linear matrix leads to a one-way nonlinearity enhancement. This is inaccessible in isotropic media which results in omnidirectional nonlinearity reduction. Thanks to the shape factor g, we can achieve both nonlinearity enhancement and directionality in this case. Another condition we study above is those linear inclusions embedded in a nonlinear matrix. It also benefited from the shape factor g that the degrees of enhancement are distinguished in different directions. When g = 0.5, our results accord with the circular particle dispersing in isotropic media. The proposed theoretical models and simulation methods can also be extended to higher-order nonlinearity, which may be expected to design flexible thermal multistability or higher heat harmonic wave generation. Taking advantage of anisotropy in configuration, directional nonlinearity enhancement can be utilized for constructing multifunctional nonlinear metadevices.

We should point out that we discuss weak nonlinearity cases in this work so that Taylor expansion can be employed, which is common in naturally-occuring solid crystals. We take first-order nonlinear terms and ignore high-order terms to clearly demonstrate the proposed methods and effects. The Taylor expansion method is universal for arbitrary order nonlinear terms, and high-order terms may be taken into consideration in some wave-like heat transfer cases. It is noted that the gain coefficient we discussed in this work has an upper limit. It depends on the intrinsic composite structures. Elliptical particles diffusing in a matrix naturally introduce this confinement of about 2.8. However, if we consider other models (for example, core-shell or diamond-shaped structure) [33, 37], the upper limits will be different. Another limitation that should be pointed out is that we employ M&G theory for deriving nonlinear effective thermal conductivities, It fits well with actual situations at the dilution limit. However, if the area fraction is large enough, simulation results will deviate from theoretical predictions. We suggest that the Rayleigh method applies to correcting the deviation between practical issues and physical models [31].

In summary, we propose an anisotropic thermal composite model for realizing directional nonlinearity enhancement. On basis of the regulation with geometrical configuration, coefficients of nonlinear terms in thermal conductivities can be enhanced in the expected directions, compared with isotropic constituent materials. By directly executing Taylor expansion on effective nonlinear thermal conductivities, we give analytical forms of nonlinearity enhancement, which is related to shape factors, linear conductivities, and area ratio of constituents. Numerical results echo the theoretical prediction and indicate the conditions for achieving directional enhancement. Moreover, we point out some limitations of our models and suggest several measures for promoting the level of directionality and enhancement. Our results may not only provide a theoretical framework for constructing directional thermal nonlinearity enhancement but also enlighten multifunctional or Janus thermal metadevice design.