In a general triplet superconductor, the spin structure of the superconducting order parameter is constrained by the underlying crystal symmetries. The structure of the UTe₂ crystals corresponds to the orthorhombic point group D _2h , and the symmetry of the odd-parity pairing states can be deduced according to the irreducible representations of the D _2h point group [30–32]. Table 1 shows the odd-parity triplet superconducting states considered in this article. Here, we do not consider the even-parity states for the D _2h point group, and we further assume strong SOC due to heavy atoms like U and Te. Due to strong SOC, orbital and spin degrees of freedom do not transform independently; instead, each space group rotation involves a rotation in spin and spatial spaces of the order parameter. In the weak SOC limit, the odd-parity states for the D _2h point group come with line nodes that are not consistent with experimental measurements.

The d-vector is real for the superconducting states listed in Table 1; hence, these states preserve TRS. These states are unitary triplet states, i.e., Δ ̂ † Δ ̂ ∝ 1 ̂ . We denote a 2 × 2 matrix in the spin space with and a 4 × 4 matrix in the Nambu-spin space with . A TRSB state is not possible with a single-component order parameter, noting that the D _2h group has only one-dimensional irreducible representations. We construct the TRSB superconducting state with a combination of two irreducible representations, and all such possible combinations are shown in Table 2. In principle, a combination of two different irreducible representations involves six real constants; however, we introduce a single-parameter model for the TRSB states. The effective d-vector is ( d 1 + i r d 2 ) / 1 + r 2 , where r is the mixing parameter that determines the relative strength of the individual order parameter. The individual d-vectors are listed in Table 1 with all the coefficients p _i=1,2..3 set to unity.

In the Nambu-spin basis, the mean field Hamiltonian reads H ̌ = ξ k σ 0 i Δ 0 d ⋅ σ σ y − i Δ 0 σ y d * ⋅ σ − ξ k σ 0 , (1) where ξ _k is the electronic dispersion, Δ₀ is the superconducting gap energy scale, and σ ₀ is the identity matrix in the spin space. We adopt a model where the electronic dispersion reads ξ k = ℏ 2 k x 2 2 m a + ℏ 2 k y 2 2 m b − μ − 2 t ⊥ ⁡ cos k z , (2) where m _a/b is the effective masses along x ̂ / y ̂ directions, t _⊥ is the hopping energy that controls the z ̂ velocity, and μ is the chemical potential. We further assume that t _⊥≪ μ. An alternate dispersion with a closed Fermi surface is considered in the Supplementary Material. The bare Green’s function is G ̌ 0 = 1 ̌ ω − H ̌ − 1 = G ̂ 11 G ̂ 12 G ̂ 21 G ̂ 22 . (3) Here, ω is the quasiparticle energy. The Matsubara Green’s function can be obtained by ω → iω _n . The 2 × 2 matrices in the spin space are G ̂ 11 = ω + ξ D ω 2 − ξ 2 − Δ 0 2 | d | 2 σ 0 + Δ 0 2 q ⋅ σ , (4) G ̂ 12 = ω 2 − ξ 2 − Δ 0 2 | d | 2 σ 0 + Δ 0 2 q ⋅ σ i Δ 0 d ⋅ σ σ y D , (5) G ̂ 21 = − ω 2 − ξ 2 − Δ 0 2 | d | 2 σ 0 + Δ 0 2 q ⋅ σ T Δ 0 i σ y d * ⋅ σ D , (6) G ̂ 22 = ω − ξ D ω 2 − ξ 2 − Δ 0 2 | d | 2 σ 0 + Δ 0 2 q ⋅ σ T . (7) Here, q = i (d ×d*), and the denominator D is D = ω 2 − ξ 2 − Δ 0 2 | d | 2 2 − Δ 0 4 | q | 2 = ξ 2 − ω 2 + Δ + 2 ξ 2 − ω 2 + Δ − 2 . (8) Here, we introduce Δ ± 2 = Δ 0 2 ( | d | 2 ± | q | ) . For the unitary case, q = 0; therefore, there is only a single energy scale. In contrast, for the TRSB nonunitary states, q ≠ 0, which leads to non-degenerate excitation energies. For single-component order parameters, q vanishes. However, for a mixture of multiple irreducible representations, q remains finite, and it can be interpreted as the spin moment of the Cooper pairs. The average of q over the Fermi surface may or may not vanish. The nonunitary states can therefore be further divided into anti-ferromagnetic (AF) and ferromagnetic (FM) states, where the average of q vanishes over the Fermi surface for the former and remains finite for the latter [9]. The six possible nonunitary states are shown in Table 2, with the possibility of nodes on a spherical or a cylindrical Fermi surface open along the z ̂ axis. For a cylindrical Fermi surface, we adopt cylindrical coordinates with k _z dependence of the gap functions replaced with sin (k _z d/2), where d is the z-axis. The factor of 1/2 is added to ensure only a single pair of point nodes in the first Brillouin zone in the unitary limit. However, this does not have any qualitative effect on our results.

In order to calculate the thermal conductivity, we need to include the effect of impurity scattering that dominates all other relaxation mechanisms at low temperatures. We consider elastic impurity scattering due to pointlike defects and include its effect through a disorder-averaged self-energy. The impurity self-energy is calculated within the self-consistent T-matrix approximation. The momentum-integrated Green’s function g ̂ has both vector and scalar components for the normal Green’s functions: g ̂ 11 = π N 0 g 0 + g ⋅ σ , (9) g ̂ 22 = π N 0 g 0 + g ⋅ σ T , (10) and the anomalous Green’s functions g ̂ 12 and g ̂ 21 vanish because the odd-parity order parameter averages to 0. In Eqs 9, 10, g is directly related to the Fermi surface average of the spin moment q and remains finite for the chiral nonunitary states, only. Using these integrated Green’s functions, we can write the T-matrix self-energy for the non-magnetic impurities as Σ ̌ = τ ̌ 3 ⋅ Σ 3 + Σ 3 ⋅ σ 0 0 Σ 3 + Σ 3 ⋅ σ T + τ ̌ 0 ⋅ Σ 0 + Σ ⋅ σ 0 0 Σ 0 + Σ ⋅ σ T , (11) Σ 3 = Γ u cot δ s cot 2 δ s − g 0 2 + g ⋅ g D imp , (12) Σ 3 = Γ u 2 ⁡ cot δ s g 0 g D imp , (13) Σ 0 = Γ u g 0 cot 2 δ s − g 0 2 − g ⋅ g D imp , (14) Σ = Γ u g cot 2 δ s + g 0 2 − g ⋅ g D imp , (15) D imp = cot 4 δ s − 2 ⁡ cot 2 δ s g 0 2 + g ⋅ g + g 0 2 − g ⋅ g 2 , (16) where δ _s ≡ tan⁻¹ (πN ₀ V _imp ) is the s-wave scattering phase shift; Γ _u = n _imp /(πN ₀); and n _imp and V _imp are the impurity concentration and impurity potential, respectively. Here, τ ̌ 3 is the Pauli matrix in the Nambu space. The τ ₃ component of the self-energy that renormalizes the electronic dispersion is ignored. It can be absorbed in the chemical potential. The impurity-dressed Green’s function reads G ̌ − 1 = G ̌ 0 − 1 − Σ ̌ = ω ̃ − ξ σ 0 − Σ ⋅ σ − Δ ̂ − Δ ̂ † ω ̃ + ξ σ 0 − Σ ⋅ σ T , G ̌ = G ̂ 11 G ̂ 12 G ̂ 21 G ̂ 22 . (17) Here, impurity-renormalized ω ̃ = ω − Σ 0 , which is obtained self-consistently. Unlike unitary superconductors, the impurity-dressed nonunitary Green’s function acquires a different structure from that of the bare Green’s function, in particular the structure in spin space for the normal component. The individual components of the Green’s function G ̌ are G ̂ 11 = L 0 + L 1 ⋅ σ D ̃ , (18) G ̂ 22 = L 0 ξ → − ξ + L 1 ξ → − ξ ⋅ σ T D ̃ , (19) G ̂ 12 = [ 2 ξ Σ ⋅ d + ω ̃ 2 − ξ 2 − Δ 0 2 | d | 2 + Σ ⋅ Σ d ⋅ σ + i Δ 0 2 q × d ⋅ σ − 2 Σ ⋅ d Σ ⋅ σ + 2 i ω ̃ Σ × d ⋅ σ ] i σ y Δ 0 D ̃ , (20) G ̂ 21 = [ − 2 ξ Σ ⋅ d * + ω ̃ 2 − ξ 2 − Δ 0 2 | d | 2 + Σ ⋅ Σ d * ⋅ σ T − i Δ 0 2 q × d * ⋅ σ T − 2 Σ ⋅ d * Σ ⋅ σ T − 2 i ω ̃ Σ × d * ⋅ σ T ] i σ y Δ 0 D ̃ , (21) where L 0 = X 0 a 0 − Δ 0 2 b 0 | d | 2 − Δ 0 2 q ⋅ Σ , (22) L 1 = X 0 + Δ 0 2 | d | 2 Σ + Δ 0 2 b 0 q − Δ 0 2 Σ ⋅ d d * − Δ 0 2 Σ ⋅ d * d . (23) Here, a 0 / b 0 = ω ̃ ∓ ξ , X 0 = b 0 2 − Σ ⋅ Σ , and the denominator D ̃ = ( ξ 2 + Q ̃ + 2 ) ( ξ 2 + Q ̃ − 2 ) , where Q ̃ ± is Q ̃ ± 2 = Δ 0 2 | d | 2 − ω ̃ 2 − Σ ⋅ Σ ± Δ 0 4 q ⋅ q + 4 ω ̃ 2 Σ ⋅ Σ + 4 Δ 0 2 ω ̃ q ⋅ Σ − 4 Δ 0 2 Σ ⋅ d Σ ⋅ d * . (24) For the unitary states, the nodes are symmetry-imposed, and for the nonunitary states, nodes may shift away from the high-symmetry directions, and their positions remain protected against disorder as long as the Σ component of the impurity self-energy vanishes. Σ can be interpreted as impurity-induced magnetization. For the chiral states, this term remains finite and gives rise to non-degenerate quasiparticle spin density and, in principle, can change the nodal positions. Here, the nodes do not refer to the zeros of the gap or the order parameter; instead, they are the zeros in the quasiparticle spectrum on the Fermi surface. In the unitary states, the gap nodes and the spectral nodes are same, unlike the nonunitary states. There are some additional triplet terms Σ ⋅ σ and (Σ ×d) ⋅ σ in Eq. 20, which reflect the impurity-induced modification of the spin structure of the Cooper pairs. It is worth mentioning that there is an impurity-induced odd-frequency pairing for the chiral nonunitary states, which is spin singlet and odd parity in nature.

After obtaining the impurity-dressed Green’s function, we calculate the electronic thermal conductivity κ using the Kubo formula that relates the thermal conductivity to the heat–current response [24]. We ignore the vertex corrections and restrict ourselves to the bare thermal–current response function. The vertex corrections are small in the strong scattering limit that is focused on in this article [33]. The diagonal thermal conductivity for a general triplet superconductor reads κ i i T = ∫ − ∞ ∞ d ω ω 2 T 2 − d n F ω d ω ⟨ N 0 v F i 2 − c 1 c 4 + c 2 c 3 + c 3 b 1 + c 1 b 2 + b 3 − c 3 + c 1 c 2 / c 4 c 3 2 + c 1 2 c 4 − c 1 c 2 c 3 ⟩ F S . (25) Here, c 1 = − 2 R e ( Q ̃ + + Q ̃ − ) , c 2 = | Q ̃ + | 2 + | Q ̃ − | 2 + 4 R e ( Q ̃ + ) R e ( Q ̃ − ) , c 3 = − 2 | Q ̃ + | 2 R e ( Q ̃ − ) − 2 | Q ̃ − | 2 R e ( Q ̃ + ) , and c 4 = | Q ̃ + | 2 | Q ̃ − | 2 , and the coefficients b _i=1,.,3 are b 1 = | ω ̃ | 2 − Δ 0 2 | d | 2 + | Σ | 2 + R e Q ̃ + 2 + Q ̃ − 2 , (26) b 2 = 1 4 | Q ̃ + 2 + Q ̃ − 2 | 2 + | ω ̃ | 2 + | Σ | 2 − Δ 0 2 | d | 2 R e Q ̃ + 2 + Q ̃ − 2 + 3 Δ 0 4 q ⋅ q + 4 | Σ | 2 | ω ̃ | 2 + R e ω ̃ 2 − 4 Δ 0 2 | Σ ⋅ d | 2 + | Σ ⋅ d * | 2 + R e Σ ⋅ d Σ ⋅ d * + 4 | ω ̃ | 2 R e Σ ⋅ Σ + 4 Δ 0 2 R e ω ̃ q ⋅ Σ * + 2 ω ̃ q ⋅ Σ , (27) b 3 = | ω ̃ | 2 − Δ 0 2 | d | 2 | α + | 2 + Δ 0 4 q ⋅ q − 2 Δ 0 4 q ⋅ q R e α + + Y Σ , (28) Y Σ = | ω ̃ 2 + Δ 0 2 | d | 2 − Σ ⋅ Σ | 2 − Δ 0 4 q ⋅ q | Σ | 2 − 2 Δ 0 2 R e α − − α − * ω ̃ q ⋅ Σ * − 4 Δ 0 2 | ω ̃ | 2 + | Σ | 2 − Δ 0 2 | d | 2 R e ω ̃ q ⋅ Σ * + 2 Δ 0 4 R e | q ⋅ Σ | 2 − q ⋅ Σ 2 + Σ ⋅ Σ q ⋅ q + 2 Δ 0 2 R e | ω ̃ | 2 − ω ̃ 2 − | Σ | 2 − Σ ⋅ Σ | Σ ⋅ d | 2 + | Σ ⋅ d * | 2 − 4 Δ 0 2 | d | 2 | ω ̃ | 2 | Σ | 2 + 4 Δ 0 2 R e α + Σ * ⋅ d * Σ * ⋅ d − 4 Δ 0 2 R e ω ̃ 2 − Δ 0 2 | d | 2 Σ * ⋅ Σ * , (29)

where α ± = ω ̃ 2 − Δ 0 2 | d | 2 ± Σ ⋅ Σ . The derivation of thermal conductivity is provided in Supplementary Material. We calculate the full temperature dependence of thermal conductivity using the self-consistently determined superconducting gap using an effective pairing potential to obtain a single transition temperature (see Supplementary Material).

We start with the single-component states based on four irreducible representations of the D _2h point group symmetry. The basis functions for these four states are listed in Table 1, where the A _1u state remains gapped unless the coefficient of one of the basis functions is set to zero. We exclude that possibility and choose the same coefficients for all three basis functions, and this choice of coefficients is adopted for the other states as well, which is a reasonable choice for a qualitative understanding of the low-energy properties. In principle, it is also possible to generate line nodes with an appropriate choice of basis function coefficients, but those possibilities are excluded considering the recent experimental results on UTe₂. Apart from the A _1u state, the B _1u state also remains gapped because the open Fermi surface along the z ̂ -axis forbids the nodes for this state. For the B _2u and B _3u states, a pair of point nodes exists along the y ̂ -axis and x ̂ -axis, respectively. Figures 1A–D show the disorder dependence of the thermal conductivity normalized to the normal state value at the transition as a function of relative reduction in the transition temperature δt _c ≡ = 1 − T _c /T _c0, where T _c0 is the transition temperature in the clean limit. For a clean system, δt _c = 0, and δt _c reaches unity as the impurity scattering kills superconductivity. The normalized residual thermal conductivity κT _c /κ _N T in the zero-temperature limit remains 0 up to a threshold value of disorder Γ _th for all four states; this Γ _th corresponds to a threshold level of T _c suppression δ t c t h . This threshold value of disorder depends on the superconducting gap structure and the strength of the impurity potential [27]. For the gapped states A _1u and B _1u , the residual thermal conductivity remains 0 for slightly higher values of disorder compared to the other two states B _2u and B _3u as expected due to the presence of impurity-induced quasiparticle states near the nodes. The normalized thermal conductivity remains very isotropic for the A _1u and B _1u states. However, for the nodal states, the residual thermal conductivity is enhanced for thermal current along the nodal directions. This trend in anisotropy in thermal conductivity continues even at finite temperatures, as shown in the temperature evolution of κT _c /κ _N T for the four single-component states in Figures 1E–H for tan θ _s = 2 and in Figures 1I–L for tan θ _s = 5, which represents the stronger scatterers. For both impurity strengths, δt _c is 0.02. In the presence of the point nodes, the thermal conductivity shows a weak maximum at very low temperature along the nodal direction for the weak scatterers, which disappears as the scattering rate increases. This is a known behavior for the superconducting states with point nodes [28], which is not present for the stronger impurity potentials.

It is useful to examine the density of states (DOS) and the structure of the low-energy quasiparticle states before discussing the thermal transport for the nonunitary states. We first report the average density of states per spin for the nonunitary states on a cylindrical Fermi surface that is open along the z ̂ -axis. Figure 2A shows the DOS for the A _1u + irB _2u state, which is a chiral state, and Figures 2B–D show Δ₋ for this state. In the gapped phase (r < 1), this state has minima near the x ̂ -axis, and a small gap is visible in the DOS. For r = 1, a pair of point nodes appear along the x ̂ -axis, and the low-energy DOS shows ω ² behavior that is expected for linear point nodes. Here, and in subsequent discussions, a “node” refers to spectral nodes, not the gap nodes. However, in contrast to a unitary state, the low-energy quadratic behavior remains confined to a very low energy scale, as compared to the unitary B _3u state. For r > 1, the nodes move away from the x ̂ -axis, and the positions of four nodes are determined by tan ⁡ ϕ = ± r 2 − 1 , where ϕ is the polar angle on the cylindrical Fermi surface. As the mixing parameter r increases, additional pairs of nodes appear in the yz plane at sin ( k z / 2 ) = ± 1 / r 2 − 1 . In the r → ∞ limit, only two point nodes along the y ̂ -axis survive, as expected for a pure B _2u state. The low-energy DOS remains quadratic in all these cases. For 1 ≤ r < 2 , the nodes remain closer to the x ̂ -axis, and for r ≥ 2 , the nodes move closer to the y ̂ -axis. The A _1u + irB _3u state also shows similar DOS to the A _1u + irB _2u state, but it has a different nodal structure. It has gap minima in the quasiparticle spectrum near the y ̂ -axis, and nodes appear along the y ̂ -axis. For r > 1, a set of four nodes appear near the y ̂ -axis at ϕ = ± cot − 1 r 2 − 1 and move toward the x ̂ -axis in the xy-plane as the value of r increases. For r > 2 , four more nodes appear at k z = ± 2 ⁡ sin − 1 ( 1 / r 2 − 1 ) in the xz-plane. Both these states are chiral and show finite quasiparticle spin density along y ̂ and x ̂ directions in the spin space. The last AF nonunitary state is A _1u + irB _1u , which is not a chiral state. It has nodes along the x ̂ and y ̂ directions for r = 1, a set of four point nodes in the xz-plane and another set of four point nodes in the yz-plane, where the k _z for the nodal position is determined by sin ( k z / 2 ) = ± r 2 − 1 . The DOS shows quadratic behavior at low energies (see Supplementary Material).

Next, we consider the ferromagnetic nonunitary states on the cylindrical Fermi surface, which are chiral states with finite Cooper pair spin moment. Figure 2D shows the DOS for the B _1u + irB _2u state along with the Δ₋ in Figures 2E–F. This state is gapped for r < 1 / 2 , and a set of four point nodes appear in the yz-plane for 1 / 2 ≤ r < 1 , whose positions are determined by sin ( k z / 2 ) = r / 1 − r 2 . For r > 1, a set of four point nodes appear in the xy-plane at ϕ = ± tan − 1 r 2 − 1 close to the x ̂ -axis and move toward the y ̂ -axis as r increases. These states show quadratic DOS at low energies. For r = 1, this state shows point nodes along the x ̂ -axis, but these are the quadratic point nodes. The first derivative of Δ₋ vanishes at the nodes for quadratic or second-order point nodes. This leads to linear DOS at low energies, as shown in Figure 2D. Similarly, for B _1u + irB _3u , a twin pair of quadratic point nodes appear along the y ̂ -axis and shows linear DOS at low energies for r = 1. For 1 / 2 ≤ r < 1 , four point nodes appear at sin ( k z / 2 ) = ± 1 − r 2 / r in the yz-plane, and for r > 1, a set of four point nodes appear in the xy-plane at cot ⁡ ϕ = ± r 2 − 1 , closer to the y ̂ -axis for r ≥ 1, which move toward the x ̂ -axis for r ≫ 1. The DOS remains quadratic, as expected. For the B _2u + irB _3u state, four point nodes are either located in the yz-plane at sin ( k z / 2 ) = ± r / 1 − r 2 for r ≤ 1 / 2 or in the xz-plane at sin ( k z / 2 ) = ± 1 / r 2 − 1 for r ≥ 2 . These states show ω ² behavior in the low-energy DOS. For 1 / 2 < r < 2 , a gap exists in the quasiparticle spectrum due to lack of nodes (see Supplementary Material). In contrast to the AF nonunitary states on the cylindrical Fermi surface, the FM nonunitary states can have at most four nodes and are expected to be more anisotropic.

We now discuss the zero-temperature limit of the thermal conductivity, which is very sensitive to the gap structure. We compare the normalized κT _c /Tκ _N along three principal directions, where the thermal conductivity along a particular direction is normalized to its normal state value at T _c along that direction. This suppresses the intrinsic anisotropy present in the electronic structure and accentuates the effect of order parameter anisotropy. Figures 3A–D show κT _c /Tκ _N in the zero-temperature limit for the A _1u + irB _1u state as a function of the relative reduction in the transition temperature δT _c . For a weakly disordered system, κ/T vanishes in the zero-temperature limit, but as the disorder level increases and crosses Γ _th , κ/T| _T→0 becomes finite and reaches the normal-state value as the superconductivity vanishes. For the A _1u + irB _1u state, κT _c /κ _N T shows isotropic behavior in the xy-plane. Note that for this state, the nodes are always along the x ̂ and y ̂ axes at the same k _z value. In contrast, the B _2u + irB _3u state shows relatively weaker level anisotropy than a unitary state. Figures 3E–G show κT _c /κ _N T for the B _2u + irB _3u state, while Figure 3H shows κT _c /κ _N T for the unitary B _3u state with nodes along the x ̂ -axis. As shown in Figure 3E, with nodes in the yz-plane κ _yy and κ _zz exceeding κ _xx , and this behavior reverses as the nodes move to the xz-plane for the B _2u + irB _3u state, as depicted in Figure 3H. As the value of r increases, the anisotropy also reduces, and for r = 1, the superconducting state becomes fully isotropic, as shown in Figure 3F. This is a special case, which has a four-fold symmetric Δ₋, leading to a very isotropic normalized thermal conductivity along the three principal directions. As r becomes larger than unity, the κ _xx starts to dominate.

Next, we look at the zero-temperature limit thermal conductivity for the A _1u + irB _2u state for two different impurity potential strengths. Figures 4A–D show κT _c /κ _N T for s-wave scattering phase shift tan⁻¹ (2), and panels Figures 4E–H show κT _c /κ _N T for θ _s = tan⁻¹ (5). κ/T| _T→0 becomes finite above a threshold disorder level, as in earlier cases. This threshold scattering rate is smaller for the stronger impurity potentials. This state has a minimum and maximum along the xz-plane, the energy gap is small in the yz-plane, and a weak maximum exists along the y ̂ -axis in the xy-plane. This excitation energy spectrum is reflected in the thermal conductivity. As κ/T becomes finite, the in-plane anisotropy is very weak with a slightly larger value along the x ̂ -axis as long as there is no node in the yz-plane. As the impurity-induced quasiparticles overcome the minima along the y ̂ -axis, the thermal conductivity along the y ̂ -axis starts to dominate. In the case of eight nodes, with four in the xy-plane and another four in the yz-plane, the thermal conductivity is always larger along the y ̂ -axis. The z ̂ -axis thermal conductivity remains close to the x ̂ -axis thermal conductivity. For the A _1u + irB _3u state, the thermal conductivity along the x ̂ -axis and y ̂ -axis shows the same behavior as κ _yy /T and κ _xx /T in the A _1u + irB _2u state, respectively.

For the B _1u + irB _3u state, in the gapped phase, i.e., r < 1 / 2 , the spectral gap on the Fermi surface is small in the yz-plane and near the y ̂ -axis in the xy-plane, which leads to larger κ/T| _T→0 along the y ̂ -axis, followed by the z ̂ direction, as shown in Figures 4I–M, for tan θ _s = 2 and tan θ _s = 5, respectively. As a set of four nodes appear in the yz-plane, the relative anisotropy remains the same, as shown in Figures 4J–N. This state also shows an elusive quadratic node for r = 1 along the y ̂ -axis. Unlike the other states with linear point nodes, for this case, κ/T| _T→0 term remains finite along the nodes, and for other directions, residual thermal conductivity remains 0 below the threshold disorder level. This state shows linear DOS at low energies, like superconductors with line nodes. Finally, the nodes appear in the xy-plane for this state as r exceeds unity and remains closer to the y ̂ -axis, and thermal conductivity along the y ̂ -axis becomes dominant, while the other two directions show very similar κT _c /κ _N T. As r increases, the in-plane anisotropy decreases, and the x ̂ -axis κ increases and becomes stronger along the y ̂ direction in the r ≫ 1 limit. For the B _1u + irB _2u state, the in-plane anisotropy found for the B _1u + irB _3u state gets interchanged.

Now, we look at the temperature evolution of the normalized thermal conductivity for the nonunitary states. At very low temperatures, the elastic scattering by the impurities is the main mechanism of relaxation. However, as the temperature increases, the inelastic scattering also becomes important, which we discuss in the subsequent section. Apart from the electron contribution to thermal conductivity, phonon thermal conductivity can also become significant. Here, we focus on electronic thermal conductivity and the effect of impurity scattering on it and the effect of underlying spectral nodes on the anisotropy in the thermal conductivity. As shown in Figure 1, the thermal conductivity shows a weak maximum as a function of temperature at very low temperatures along the nodal directions. There is no evidence for such a feature in the experimental measurements [7, 34]; therefore, we set θ _s = tan⁻¹ (5) for the rest of our discussion (see Supplementary Material for the weaker impurity scatterers). We first consider the A _1u + irB _1u state, which shows very isotropic residual thermal conductivity, and the thermal conductivity remains isotropic as a function of temperature, as shown in Figures 5A–D. For the clean system with δt _c = 0.02, at very low temperatures, the thermal conductivity remains negligibly small, and as the temperature increases, κ/T increases T ² at very low temperatures. This behavior is observed for both nodal and gapped systems. Note that in superconductors with point nodes, there are very few states available at the Fermi energy, as shown in Figure 2. For dirtier systems such as δt _c = 0.15, there are sufficient quasiparticle states at the Fermi level to provide nonzero thermal conductivity, which is depicted in Figures 5E–H. At low temperatures, κ/T remains independent of temperature, and as it increases, once more, quasiparticles become relevant for transport as the temperature increases.

Next, we consider the A _1u + irB _2u state, which shows enhanced thermal conductivity along the y ̂ -axis. The zero-temperature trends in anisotropy continue as the temperature increases, as shown in Figure 6. For δt _c = 0.02, κ/T is close to 0 and increases as temperature increases beyond a critical value. For the gapped case (r < 1) and for the x ̂ -axis nodes (r = 1), κ/T is slightly larger than that in the other two directions, but as the temperature increases, the other two directions begin to increase and dominate because the temperature overcomes the gap minima in the yz-plane. For r > 1, there are four point nodes that move toward the y ̂ axis, and it is reflected as a larger κ/T along the y ̂ -axis. For r > 2 , four more nodes in the yz-plane make κ/T along y ̂ and z ̂ directions larger than that in the x ̂ direction. For dirtier systems, thermal conductivity becomes finite and T-independent at very low temperatures, as in the earlier case, but the anisotropy remains similar to its zero-temperature limit. For A _1u + iB _3u , the behavior of x ̂ and y ̂ directions gets interchanged (see Supplementary Material).

Next, we consider the B _2u + irB _3u state, which is one of the FM nonunitary chiral states. For this state, there are four point nodes in the yz-plane for r < 1 / 2 , and for r > 2 , there is a set of four point nodes in the xz-plane. Figure 7 shows the thermal conductivity as a function of temperature, and as expected, κ/T along the nodal directions dominates. For r < 1, the thermal conductivity is enhanced along the y ̂ and z ̂ directions, and as r increases, the system becomes gapped with enhanced thermal conductivity in the y ̂ and z ̂ directions, but the anisotropy reduces. The state r = 1 shows completely isotropic behavior, and as r exceeds unity, this state shows a x ̂ -axis-dominated thermal response, as the nodes reappear in the xz-plane. The anisotropy remains qualitatively same as the temperature increases. Finally, we consider the B _1u + irB _3u state, which is another possible FM nonunitary state that shows enhanced thermal conductivity along the y ̂ direction, as shown in Figure 8. Here, the anisotropy changes significantly as the temperature increases. At low temperatures, the y ̂ -axis dominates due to its vicinity to the point nodes; however, as the temperature increases, κT _c /κ _N T increases along the x ̂ directions. This happens because at lower temperatures, the lower-energy branch of the quasiparticle excitations ξ k 2 + Δ − 2 dominates, which has more quasiparticle states along the y ̂ directions, but at higher temperatures, the ξ k 2 + Δ + 2 branch of the quasiparticle excitations begins to contribute, which, for this state, has minima along the x ̂ -axis. The overall anisotropy for this state is relatively less compared to other states, except those that show fully isotropic behavior as a function of temperature or disorder. For the B _1u + irB _2u state, κ/T along the x ̂ and y ̂ directions gets interchanged by the anisotropy shown by the B _1u + irB _3u state.

As mentioned in the previous section, at very low temperatures, i.e., T ≪ T _c , the elastic scattering from impurities is the only mechanism that determines the scattering rate. However, at higher temperature, inelastic scattering from a bosonic mode is possible. We consider a simple scenario where there is a dispersionless bosonic mode that couples with the fermions with an effective coupling constant g _fb . We further assume that the coupling does not depend on the spin degree of freedom of the fermion. The lowest-order self-energy for the fermions reads Σ i n i ω n , k = g f b 2 T ∑ m , q G i ω n − i Ω m , k − q D i Ω m , q . (30) Here, D is the bosonic Green’s function and ω _n and Ω _m are the fermionic and bosonic Matsubara frequencies, respectively. After performing the Matsubara summation, we obtain Σ i n ω = − g f b 2 2 π ∫ − ∞ ∞ d x ∫ − ∞ ∞ d y N y D ′ ′ x x + y − ω − i η coth β x / 2 + tanh β y / 2 . (31) Here, the real part of the self-energy contributes to mass renormalization, and the imaginary part modifies the scattering rate, which is a function of quasiparticle energy and temperature. Noting that UTe₂ has very high effective mass, we ignore the real part of the self-energy. The imaginary part of the self-energy is Σ i n ′ ′ ω = − g f b 2 2 ∫ − ∞ ∞ d x N ω − x D ′ ′ x coth β x / 2 + tanh β ω − x / 2 , (32) where D ′ ′ ( x ) ≡ x / ( x 2 + Ω 0 2 ) is the bosonic DOS and Ω₀ is the characteristic energy scale associated with the bosonic mode. In the context of UTe₂, we expect Ω₀ ≫ T _c ; therefore, we approximate the bosonic DOS as D ′ ′ ( x ) ≈ x / Ω 0 2 . In the zero-temperature limit, Σ i n ′ ′ ∝ ω n + 2 , where the DOS for the superconducting state behaves like ω ⁿ at low energies. For linear point nodes, Σ i n ′ ′ ∝ ω 4 , and for quadratic point nodes or line nodes, it behaves like ω ³. Similarly, in the static limit (ω → 0), the imaginary part of the inelastic self-energy reads Σ i n ′ ′ ω = 0 , T = − 2 g f b 2 ∫ 0 ∞ d x N x D ′ ′ x sinh ⁡ β x . (33) Since the thermal response integrand is peaked at ω = 0, we only retain the temperature dependence of the inelastic self-energy at ω = 0. At very low temperatures, the inelastic scattering rate behaves likes T ⁴ for the point nodes and T ³ for the line nodes in the static limit. This is sufficient to understand the qualitative effect of inelastic scattering. The prefactor in Eq. 33 is fixed by the value of the inelastic scattering rate at T _c .

Figure 9 shows the effect of inelastic scattering on the electronic thermal conductivity. We show the results for two cases as the rest of the cases are qualitatively similar. The first row of Figures 9A–D shows the normalized thermal conductivity for the A _1u + irB _2u state, which shows enhanced thermal conductivity along the y ̂ -axis. Note that the normal state thermal conductivity includes the inelastic scattering. We fixed the value of r at 1.2, which leads to four point nodes in the xy-plane. The primary effect of inelastic scattering is the formation of a peak below T _c . This peak appears first for the directions that have larger thermal conductivity, as shown in Figure 9A for δ = 0.02. As the system becomes dirtier, the peaks get smeared, as depicted in Figure 9B for δt _c = 0.15. If the strength of inelastic scattering increases, the peaks also strengthen, as shown in Figure 9C, and it may survive in dirtier systems, as illustrated in Figure 9D. Similar trends continue for the other states as well. Figures 9E–H show the normalized thermal conductivity for the B _1u + irB _3u state, which shows the same qualitative behavior as a function of disorder and inelastic scattering strength. However, this state shows a nonmonotonic variation of anisotropy, and the peaks for different directions remain quite close to each other.