Influence of ferroelastic domain walls on thermal conductivity

Enabling on-demand control of heat flow is key for the development of next-generation electronic devices, solid-state heat pumps, and thermal logic. However, precise and agile tuning of the relevant microscopic material parameters for adjusting thermal conductivities remains elusive. Here, we study several single crystals of lanthanum aluminate (LaAlO$_{3}$) with different domain structures and show that ferroelastic domain walls behave as boundaries that act like efficient controllers to govern thermal conductivity. At low temperature (3 K), we demonstrate a fivefold reduction in thermal conductivity induced by domain walls orthogonal to the heat flow and a twofold reduction when they are parallel to the heat flow. Atomistic calculations fully support this experimental observation. By breaking down phonon scattering mechanisms, we also analyze the temperature dependence of the thermal conductivity to derive a quantitative relation between thermal conductivity variations and domain wall organization and density.

Their efficiency is parametrized through the ratio (R) of their high (κhigh) and low (κlow) thermal conductivity states.Some are based on fluids, but solid-state thermal switches operating through conduction mechanisms are more appealing because of their resilience and compactness.For example, a single material such as VO2 exhibits a switching ratio R~1. 6, even though it works in a temperature range limited by the temperature of its transition between thermally insulating and conductive states [9,10].Phase-change materials exhibit larger thermal conductivity differences depending on their crystallinity [11], but lack a demonstration of reversible switching.Alternatively, multilayers of ferromagnetic and non-magnetic conductive layers can exhibit thermal conductivity ratios R~1.8 under a magnetic field [12].However, the main limitation of solid-state thermal switches operating through conduction mechanisms remains their comparatively low thermal conductivity ratios (R < 1.8).
Recently, ferroic oxides have been proposed for obtaining fast and reversible solid-state thermal switches with large switching ratios in broad temperature ranges [13][14][15][16].They rely on interactions between phonons conducting heat, and spontaneously occurring topological defects known as domain walls [17] that move in response to an electric field [18] or uniaxial stress [19].Molecular dynamics simulations reveal that the thermal conductivity could be divided by up to R=3.7 under the action of ferroelastic domain walls [14], and first-principles calculations show a thermal conductivity ratio R=2 when the density of domain walls is increased with an electric field [15].Seminal experimental works, where the precise orientation and density of domain walls are unknown, mention the influence of domain walls on the thermal conductivity at low temperature [20][21][22].More recently, at room temperature, thermal conductivity variations as high as R=3 were demonstrated in ferroelectric thin films with different densities of domain walls [23], and switching ratios as high as R=1.2 where obtained under the application of an electric field [24][25][26].
Nevertheless, some published results question the influence of domain walls on the thermal conductivity, e.g. in bismuth ferrite [27,28].A reason behind these discrepancies is that most measurements are performed only on thin films at room temperature, with complex domain structures, where the influence of residual strain, substrate and defects can be substantial.Furthermore, the thermal conductivity is assessed by thermoreflectance, after removing the response from the interface with the substrate that dominates the initial signal, and only in one direction, missing out anisotropies induced by different orientations of domain walls.These anisotropies are also excluded from current models relying on an interfacial thermal resistance at the wall [23,27].As such, questions of how domain walls influence the thermal conductivity and how their organization and density should be tuned to have a large impact are still fully open.
In this paper, we investigate the thermal conductivity of bulk single crystals of ferroelastic lanthanum aluminate (LaAlO3) with different density and orientation of domain walls.Here, we measure the thermal conductivity of five samples in different directions and in a broad temperature range -from room temperature down to 2 K.In this temperature range, i.e. far from the phase transition occurring near 850 K, domain walls in LaAlO3 have an estimated thickness of 20 Å [29] and are pinned [30,31], which means that the domain structure is stable and does not change.This unique set of data reveals large thermal conductivity variations below 10 K when the density of domain walls orthogonal to the heat flow is varied, leading to a fivefold reduction in thermal conductivity.It also reveals that domain walls have an influence on heat flows even if they are parallel to them.Still, the thermal conductivity of the material depends of the direction of the heat flow with respect to domain walls, with a twofold difference in thermal conductivity between both directions.
We account for these results by analysing the full temperature dependence of the measured thermal conductivity within the framework of the Holland model by considering all relevant phonon scattering mechanisms.We are thus able to relate quantitatively the evolution of the thermal conductivity to ferroelastic domain walls densities and orientations, successfully described by the Casimir limit, here calculated numerically to account for finite size effects.The influence of domain walls on the thermal conductivity is also computed using atomistic simulations of the heat transport (both parallel and orthogonal to the walls) based on second-principles models.By demonstrating that heat flows can be effectively controlled by domain walls, these results provide a general approach to model the influence of these planar defects on thermal conductivity, and pave the way for designing efficient, microscopic thermal switches.

A. Ferroelastic domains engineering
Figure 1 shows images of the 5 commercial single crystals of LaAlO3 studied, whose surfaces are orthogonal to the [001]pc direction ("pc" stands for pseudo-cubic), obtained with an optical microscope working in reflection (images in transmission in Supplementary Note 1 [32]).Sample F [Fig. 1(a)] has been annealed in air at 1200 K for 6 h with a heating and cooling rate of ~10 K min -1   to obtain a monodomain state, i.e. free of domain walls, which can be achieved thanks to a redistribution of defects at high temperature that act as pinning centres for domain walls [33].Since there are no domain walls, the sample appears as uniform both in reflection and transmission [Fig.1(a) and Fig. S1].Note that faint structures optically visible are not ferroelastic domains but negative traces of previous domain patterns that remain as a surface relief [34].Fig. 1(b) (sample A) and 1(c) (sample C) are as-received single crystals with a very regular domain pattern, clearly visible because of differences in reflectivity between ferroelastic domains.Domain walls are strictly orthogonal to the surface, indicating that they are planes orthogonal to the [010]pc direction, in agreement with the literature [35,36].Domain structures observed in Fig. 1(e) (sample B) and Fig. 1(f) (sample D) result from quenches of single crystals from 680 K to a room temperature silicon-oil/water mixture implying a large increase in the number of domain walls, which remain mostly orthogonal to the [010]pc direction.For all samples, the direction of the ferroelastic distortion in the domains is indicated in Fig. 1(d).The size of the domains, and hence the density of domain walls, has been measured directly from optical images and is shown in Fig. 1(g) for all samples.Samples A and C exhibit analogous distributions, with a mean sizes around 19 µm +/-2 µm and 15 µm +/-2 µm, respectively.Samples B and D exhibit smaller domains, and thus higher density of domain walls, with a mean size around 9 µm +/-2 µm.The number of domains is smaller in samples A and B, whose long lengths are along the [100]pc direction, compared to samples C and D, but the volume occupied by domain walls remains identical.

B. Thermal measurements
The temperature dependence of the measured thermal conductivity κ is shown in Fig. 2(a).It is probed in directions where the ferroelastic distortion of domains is the same by symmetry [Fig.1(d)], and thus cannot influence the thermal conductivity.The three regimes of transport usually observed in insulating oxides are recovered: the low temperature boundary scattering region characterized by a power law T-dependence, the presence of a high maximum near 25 K that depends on the effects of impurities, followed by a decrease at higher temperatures due to Umklapp scattering.At room temperature, thermal conductivity values lie between 11 and 13 W m -1 K -1 , in agreement with literature [21].In the boundary scattering regime, a decrease of the thermal conductivity is observed in samples with domain walls parallel to the heat flow (A and B), which is even stronger when domain walls are orthogonal to the flow (C and D).This is illustrated in Fig. 2(b) by the strong enhancement of the switching ratio R achieving values as high as 5 near 2 K and demonstrate the effect of domain walls on thermal transport.We verified that annealing and quenching do not influence the thermal conductivity if the samples stay monodomain (Supplementary Note 3 [32]).In particular, in the optical images of quenched samples some cracks are visible (Supplementary Note 1 and 3 [32]).Even though these cracks are mostly at the surface and do not extend across the thickness of the samples, except on the edges, they are particularly visible because the images are taken in transmission.We verified that these cracks have little influence on the thermal conductivity by comparing in Fig. S2 [32] the thermal conductivity of two monodomain samples: one obtained after annealing (sample F) and one obtained after quenching.At 3 K, the thermal conductivity of the quenched sample is reduced by a factor 1.4 only, far from the fivefold reduction in thermal conductivity observed between samples F and D.
To further analyze the influence of domain walls, we performed specific heat measurements down to 2 K [Fig.2(c)].Below 10 K the usual low temperature Debye regime is recovered with a welldefined T 3 behavior, characteristic of acoustic phonons contribution.We calculate the Debye temperature TD=375 K according to the expected contribution C(T<<TD)=12 4     , is unveiled in Fig. 2(d).In the low temperature boundary scattering regime it saturates differently depending on the orientation and density of domain walls.In particular, L is more constrained when domain walls are orthogonal to the heat flow and their density is high, as in sample D compared to sample C, rather than in the case of domain walls parallel to the heat flow as in samples A and B. In the latter case, it is recovered that a higher domain-wall density in sample B reduces more efficiently the mean free path than in sample A. Whichever the orientation of domain walls, L is always lower than in the monodomain sample F.

A. Simulations of heat transport with a domain wall
We first compare our results with atomistic heat transport simulations.To this end, we derive a second-principles potential for LaAlO3 (see Appendix B) and we use it to optimize cells with domain walls along the {001}pc direction, which we find to be ~25 Å wide in the low temperature limit, close to the experimental value of 20 Å (ref.[29]).We then make use of non-equilibrium Green's function techniques to compute the thermal conductance of LaAlO3 in three configurations: monodomain state, with domain walls orthogonal and parallel to the heat flow [Figs.3(a) and 3(b)].We observe that domain walls reduce the thermal transport moderately when parallel to the heat flow and do so more markedly when orthogonal to it.This agrees very well with our experimental results.
In the simulation, the fact that the difference in thermal conductances increases with temperature indicates that the thermal resistance is attributed mostly to optical modes.This is corroborated by the simulated frequency-resolved transmission curves [Figs.3(c) and 3(d)], which show that optical modes of low and medium frequencies (3-20 THz) undergo a stronger scattering (and the scattering is more pronounced in the orthogonal configuration).In contrast, the experiments show that acoustic modes dominate at low temperature.The discrepancy may arise from the fact that the real (experimental) walls in the multidomain configuration induce a local strain, especially in the high walldensity limit, while in the simulations we considered a relatively sparse regime.This strain could induce large scattering in the acoustic modes, and hence reduce further the thermal conductivity.
Also, the experiments could be measuring partly additive scattering by consecutive walls if the scattered phonons recover some of their population through anharmonic interactions before hitting the next wall.This effect would not be captured by the simulations since they are carried out in the harmonic approximation and do not account for phonon-phonon scattering.For that same reason, the simulated thermal transport curves increase monotonically as a function of temperature, as phonons become increasingly populated.Including anharmonic scattering would recover the correct high-temperature behavior, where the effect of the domain walls vanishes because phonon-phonon scattering becomes dominant but comes at a prohibitive computational cost given the sizes of the simulation boxes employed.Still, our simulations give the right trends for the thermal conductivity reduction caused by domain walls in LaAlO3.

B. Thermal conductivity of the monodomain sample
In a more general framework, the thermal conductivity can be described by using the semi-classical Boltzmann equation in the frame of the relaxation time approximation  for phonons characterized by Bose-Einstein statistics with zero chemical potential.It is written as the integral over the acoustic phonon frequencies by using a Debye model [38]: It follows in the case of a boundary scattering, for which the relaxation time is neither frequency nor temperature dependent (=0 and  =B), that the thermal conductivity varies as T 3 and κ = For the monodomain sample of LaAlO3, this regime is never reached since the measured thermal conductivity varies as T 2 at low temperature in Fig. 2(a).Already observed in some perovskite oxides [20], this behavior is the signature of phonon scattering by dislocations [39] with a linear frequency dependence of  -1 , namely with  = -1 and no temperature dependence of the relaxation time.Therefore, it follows from Eq. ( 2) that  varies as T 2 at low temperature, which explains why the mean free path in Fig. 2(d) is not constant.Within the model of Holland [40], which considers two types of phonon polarization with 1 longitudinal (L) and 2 transverse (T) modes, κ = where κ  and κ  are both defined as in Eq. ( 1) and the factors 1/3 and 2/3 originate from the numbers of longitudinal and transverse modes, respectively.Additional scattering processes are also considered within the model by including crystalline boundaries, impurities through mass differences, the three-phonon Normal processes (L and T) as well as the Umklapp ones (L and T), and dislocations, such that the relaxation time entering Eq. ( 1) is where each j is the relaxation time for a single scattering mechanism [40].
The transverse contribution is further decomposed as κ  = κ  + κ  in order to take into account the specific frequency range of transverse Umklapp scattering [40], which has been adjusted numerically to fit the thermal conductivity.The best adjustment is reached here with a transverse cutoff frequency D/1.5 as summarized in Table 1.As shown in Fig. 4(a), the thermal conductivity measured in the monodomain sample F is successfully described by these 3 contributions, L, TN and TU.According to their definition in Table 1, the total relaxation times used for the calculations are:

C. Thermal conductivity of multidomains samples
We have refined the parameters listed by calculating numerically the integrals defined in Eq. ( 1) at each temperature and for each component to fit the measured thermal conductivity in all the samples.The found factors Ai (Table 1) agree quantitatively with the ones reported for several insulating oxides [41].One must also emphasize that the only parameters that have been varied from one sample to the other are the dislocation factor and the boundary length LB, which depends on the density and orientation of ferroelastic domains in the samples.Since dislocations act as nucleation or pinning centers for domain walls, one may assume that the contribution from the walls likely takes over the contribution from dislocations.This could then explain why it has been found that the dislocation factor decreases from 2.6 × 10 -5 (sample F) down to 1.15 × 10 -5 (samples A and C) and 0 (samples B and D) when the domain-wall density increases.While it remains difficult to ascribe a precise dislocation density from the AD factor [42], it is known to be proportional to the dislocation density [42,43], and similar values have been found in materials with dislocation densities ranging from 1 × 10 7 up to 5 × 10 7 cm -2 , which seems here a reasonable order of magnitude.

D. Specularity parameter
Boundary scattering has been described here within the framework of Casimir [44,45].It assumes that the temperature is so low that phonons interact only with the boundary and that the scattering is completely diffuse.This implies that incident phonons which are absorbed by the various surfaces of the sample are then re-emitted with the equilibrium distribution corresponding to the local temperature.However, it appears in Fig. 2(d) that the mean free paths still increase at low temperatures in contrast to the expected constant boundary regime behavior.This is because the boundary scattering of phonons as described by Casimir assumes perfectly diffusive surfaces where phonons are scattered with equal probability into any directions, which is only ensured if the roughness of the surface is higher than the phonon wavelength.Nevertheless, the latter increases when the temperature is lowered, and the previous diffusive criterion is inevitably no longer valid for a sufficiently smooth surface, as shown in Fig. 4(b).In this case, phonons undergo specular reflection on the boundary rather than diffuse scattering and their effective mean free path increases beyond the boundary length usually referred to as the Casimir limit.This effect has been described by Ziman who has proposed to account for it through an average specularity parameter p dependent on the roughness () distribution P() and the phonon wavelength  (ref.[46]): Accordingly, an effective mean free path   eff is related to the boundary length LB in Eq. ( 3) by considering multiple specular reflections, which lead to a strong enhancement if  → 1 (perfectly smooth surface) or to the bare boundary length if  = 0 (diffuse scattering).Whereas the roughness distribution is a priori unknown, one may assume that for a polished surface this is a strongly decreasing function of  which can be approximated by an exponential form such as (η) ≈  −η η 0 ⁄ η 0 (ref.[47]).This specularity parameter becomes then  ≈ 1 −  −λ 4 ⁄ πη 0 and   eff = (2 λ 4 ⁄ πη 0 − 1)  , by considering the dominant phonon wavelength λ ≈  2).

E. Casimir limits from domain distributions
To discuss quantitatively the role of the density and orientation of domain walls, the optically measured sizes of domains must be converted into a relevant boundary length that can be precisely compared to the boundary length LB (or LB') inferred from the analysis of the thermal conductivity.
For this, as explained in Appendix F, we have calculated the so-called Casimir limit which accounts for the finite size effects in the samples (e.g.size of ferroelastic domains) as a function of their length c in Fig. 4(d), by using their experimental dimensions (width a=1600 µm, thickness b=500 µm).
For small rectangular rods of extension less than 100 µm it appears that the Casimir limit is mainly given by one configuration factor (see Appendix F) such as 3Fzz/2 according to Eq. (A6), due to the smallness of the lateral surfaces.This factor will thus be used to analyze the effect of domain walls orthogonal to the heat flow.The variation of the Casimir limit as a function of the width a in the inset helps to understand the effect of domain walls parallel to the flow because in this case the limit of infinite rod appears justified.In fact, we find that each domain can be considered as an independent rectangular rod with its own thermal conductivity.If walls are orthogonal to the flow, the sample can be represented by a collection of domains in series from the transport point of view, whereas if walls are parallel to the flow, domains must be considered in parallel.It follows that in the former case thermal resistances are additive whereas in the latter thermal conductances are additive.By writing the total thermal conductivity in the boundary regime as κ ⊥,∥ = 1 3    ⊥,∥ , depending on the walls' orientation (orthogonal or parallel), with the domain sizes ci or ai respectively, the macroscopic Casimir limit  ⊥,∥ is then related to the distribution of domains through their size and respective Casimir limit LC,i: The resulting macroscopic Casimir limits LC are summarized in Table 2 and compared to the boundary lengths LB (and LB') found experimentally from the fit in Fig. 2(d).

IV. CONCLUSIONS
Thermal conductivities of the monodomain sample and the samples with domain walls orthogonal or parallel to the heat flow [Fig.2(a)] are successfully reproduced in the frame of the Holland model, by using boundary lengths indicated in the caption with a specular enhancement factor with a rms roughness =7 nm (or 4 nm for the frequency-dependent specularity).Accordingly, the mean free path of the 5 samples appears suitably described in Fig. 2(d) with boundary lengths originating from different domain sizes in agreement with the inferred Casimir limits for samples F, A and B (LB≈LC) or consistent with   ≈   × (1.7 ± 0 .1)for samples C and D. These results demonstrate that domain walls behave as boundaries that act as thermal conductivity controllers according to their orientation.The distribution of domain walls implies a distribution of thermal conductivity considered either in series or in parallel depending on the orientation of the walls with respect to the heat flow.
The fact that the found boundary length in the orthogonal configuration gets closer to the Casimir limit if one considers a frequency-dependent specularity parameter likely suggests that a more precise description of domain walls, including roughness distributions, is required in these conditions of very short domain sizes.
Overall, we demonstrated that tuning the number and orientation of ferroelastic domain walls is a powerful way to govern thermal conductivity, and we provided a general approach to model the influence of these planar defects on thermal conductivity, paving the way for efficient and compact thermal switches in ferroic materials.At room temperature, LaAlO3 exhibits a rhombohedral structure, with space group 3 ̅  (ref.[49]).

ACKNOWLEDGMENTS
It thus exhibits eight domain states: four related to the loss of point-group symmetry and four related to the loss of translational symmetry (antiphase boundaries) [34,50,51].The former leads to six possible pairs of domains.For each of the six pairs there are two possible orientations of domain walls (Supplementary Note 2 [32]).The domain structure was observed with an optical microscope (Olympus BX60) operating in reflection with a polarizer and an analyser, and a x10 objective.
Quenches of single crystals from 680 K to a room temperature silicon-oil (Fisher Scientific)/distilled water mixture (~10/90 volume ratio) came after several attempts at different temperatures with other quenching mediums.Typically, quenching in air and liquid nitrogen did not affect the domain structure, contrary to quenching in oil and water.However, the latter led to the formation of significant cracks in the samples, hence the choice of an oil/water mixture as the best option to reach high density of domain walls without breaking samples.

APPENDIX B: NUMERICAL SIMULATION OF THE THERMAL CONDUCTANCE REDUCTION IN THE HARMONIC APPROXIMATION
We derive a second-principles [52,53] model for LaAlO3 based on first-principles simulations.
The density functional theory calculations used to fit the data were carried out using the Vienna Ab initio Simulation package (VASP) [54,55], making use of the PBEsol [56] implementation of the generalized gradient approximation for the exchange-correlation functional.The model can be viewed as a multi-variable Taylor series of the potential energy of LaAlO3, taking the perovskite cubic phase as the reference structure and treating all atomic distortions explicitly; we work with an effective potential that includes harmonic and anharmonic interatomic couplings up to fourth order.Domain walls were optimized by seeding an abrupt domain wall in simulation cells elongated in the direction perpendicular to the domain wall and performing Monte Carlo simulated annealings starting from 10 K, which yield smooth profiles for the octahedra rotation pattern.Two domain walls are included in the simulation box to conform with periodic boundary conditions.
The thermal conductance was calculated using non-equilibrium Green's function techniques as described in ref. [57], where the second-order interatomic force constants (IFCs) were computed with the help of the phonopy package [58].For the configuration with domain walls orthogonal (respectively parallel) to the heat transport (taken as the [100]pc direction), a simulation box of 28x2x2 (respectively 3x14x2) pseudocubic 40-atom unit cells was used to compute the IFCs.The latter were Fourier-transformed in k space over the [010]pc and [001]pc directions, and the thermal conductance along [100]pc was evaluated with a supercell of 28x1x1 (respectively 3x14x1) over a grid of 9x9 (respectively 1x15) transverse k points.The contacts were defined at both extrema of the supercell with symmetric regions made of 2x1x1 (respectively 1x14x1) pseudocubic cells.

APPENDIX C: SPECIFIC HEAT AND THERMAL CONDUCTIVITY MEASUREMENTS
The specific heat has been determined with a calorimeter of a Physical Properties Measurements System from Quantum Design using a relaxation method with a temperature rise of 2% of the sample temperature.Each measurement has been duplicated without the sample to compensate for the temperature dependence of the specific heat of the grease used for a good thermal contact between the sample and the calorimeter platform.The parameter indicating the quality of the thermal contact between the sample and the platform, the so-called sample coupling, has remained between 98% and 100% during each measurement, ensuring thus the reliability of the results.The sample measured was a small piece of 1.7x1.7x0.5 mm cut from the same single crystal as sample B.
To perform thermal conductivity measurements, 100 nm gold electrodes have been deposited through a mask to ensure low contact resistances by using silver paste.Thermal conductivity measurements have been performed under a maintained secondary vacuum (P < 10 −5 mbar) by using a four points configuration with a Physical Properties Measurements System from Quantum Design and temperatures rise of 1%.All the measurements have been performed at decreasing and then increasing temperatures to avoid thermal hysteresis by averaging both measurements.

APPENDIX D: DEBYE MODEL OF THE SPECIFIC HEAT
We emphasize that the Debye model here used accounts for the 3 expected acoustic phonon modes, 1 longitudinal and 2 transverse, which lead (as expected from a physical point of view) to a high temperature asymptotic specific heat C(T>>TD)=3 NavkB of the order of 25 J K -1 mol -1 , the rest being ascribed to the 12 optical phonon modes, which do not contribute at low temperatures.This explains the difference between the inferred Debye temperature here and those that can be found in literature where the authors frequently use an effective temperature-dependent Debye temperature to account for the overall variation of the specific heat.Then, their values need to be divided by Here, dSxy is an element of the cross section and the sum runs over all the surfaces of the sample by involving the corresponding configuration factor   = 1  ⁄ ∬ cosθ  cosθ     2 ⁄     with rzi the distance between dSxy and dSzi, xy and zi the angles between rzi and the normal vectors to the elements dSxy and dSzi (ref.[43][44][45]).Also, it is assumed that the sound velocity v is isotropic.It follows then that the classical relation between the thermal conductivity and the mean free path in the In the case of an infinite circular cylinder, LC is equal to its diameter.If it is an infinite square rod LC is about 1.12 times the side of the square [44,45], and Harrison and Pendrys [59] have provided an analytical expression in the case of an infinite rectangular rod of thickness b and width a which is slightly simplified below, with the ratio r=b/a.
In the currently interesting case of the finite rectangular rod, the Casimir limit can be further simplified by symmetry as   = )) (A6) The other configuration factors Fzx and Fzy are symmetric, and the former can be inferred from the latter according to Eq. (A7) with the replacement of b by a and vice versa.Since the last term in the following equation appears difficult to integrate analytically in contrast to the first two terms, Fzy remains below as an integral form which can be numerically computed.
The use of Eq. (A6) and Eq.(A7) allows then to compute the Casimir limit for finite size samples.compared to the rms roughness 0=7 nm.The intersection of both curves defines a cross-over from a purely diffusive regime at high temperatures and a low temperature one where boundary reflections become specular.(c) It follows that the mean free path L increases by exceeding the boundary length when the temperature is lowered, as for the sample D with its inferred mean free path L compared to the calculated ones with and without specular reflection.The inset displays the variation of the boundary lengths including (LB') or not (LB) a frequency dependent specularity parameter p(), as a function of the calculated Casimir limit related to the 5 investigated samples.(d) Variation of the Casimir limit as a function of c according to Eqs. (A4), (A6) and (A7) for a finite rectangular rod compared to the result obtained with Eq. (A5) for an infinite one [59].The inset displays the variation of the Casimir limit   ∞ as a function of the width a according to the latter equation.

/ 5
NAv kB (T/TD)3 , with NAv the Avogadro's number and kB the Boltzmann constant.By using the mass density M/V=6.52 × 10 6 g m -3 the molar mass Mmol = 213.88g mol-1 , and  −1 , in agreement with literature[21,37] (ℏ = ℎ 2π⁄ is the Planck constant).The temperature dependence of the phonon mean free path L, calculated as κ = 1 3 Here, frequencies are made dimensionless by introducing variables  = ℏω k B  ⁄ and   = ℏω k B   ⁄ .The relaxation time depends on phonon scattering mechanisms and is in most cases a function of frequency and temperature such as  =   ( with a possibly T-dependent   and an exponent  characteristic of the considered scattering process,  being a related characteristic frequency.At low temperature (typically T < TD/20), the previous integral can be extended to infinity and in the simple case of one type of scattering, the thermal conductivity can be calculated analytically by introducing Gamma  and Riemann zeta  functions:

1 3𝐶
with the temperature independent mean free path LB = B.

ℎ𝑣 2 .
821k B  .As an example, a comparison is made in Fig.4(c) between the mean free path inferred from thermal conductivity measurements in sample D without specularity, and the expected ones for a boundary length of 15 µm with a rms roughness  = 7 nm.The effect of specular reflections appears below ~10 K which corresponds to the expected crossover in Fig.4(b) between the high temperature diffuse scattering regime (when ), and the low temperature specular regime (when >).It follows that instead of saturating up to the boundary length, the mean free path still increases if temperature decreases [Fig.4(c)].Another comparison is also made by considering a frequency-dependent specularity parameter[48] (as explained in Appendix E) to take into account the different effects of the roughness depending on the short or long phonon wavelengths.Shorter and likely more realistic boundary lengths   ′ are then inferred (Table This paper was cofunded by the European Union [European Research Council (ERC), DYNAM-HEAT, No. 101077402].Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the ERC.Neither the European Union nor the granting authority can be held responsible for them.Work at LIST was supported by the Luxembourg National Research Fund through Project No. FNR/C21/MS/15799044/FERRODYNAMICS.We acknowledge financial support by MCIN/AEI/10.13039/501100011033under grant PID2020-119777GB-I00, and the Severo Ochoa Centres of Excellence Program under grant CEX2019-000917-S, and by the Generalitat de Catalunya under grant 2021 SGR 01519.Low-temperature measurements exploring the effect of domain walls on thermal conductivity in LaAlO3 were also made concurrently at Belfast by Bogdan Zhigulin, Fran Kurnia, Marty Gregg, and Raymond McQuaid.This led to initiation of the second-principles modeling study by HA, MR, RR, and JI.APPENDIX A: DOMAIN ENGINEERING Single crystals of LaAlO3 10x10x0.5 mm were obtained from Pi-Kem Ltd, with surfaces orthogonal to the [001]pc direction and edges along {001}pc directions.They were cut with a diamond wire saw (Well 3500 Precision) in smaller pieces (~1.6x8x0.5mm),keeping edges along {001}pc directions.Their main surfaces (~1.6x8mm) were optically polished by the supplier, leading to RMS roughness below 1 nm (measured by atomic force microscopy).Samples F, A, C and D come from the same single crystal of 10x10x0.5 mm.
5 1/3 due to the 5 atoms in the unit formula LaAlO3 to recover the Debye temperature found here.Both models lead to the same value when this extra term is compensated for in the calculation of the sound velocity.The experimental determination of TD allows to compute numerically the full temperature dependence of the Debye contribution as defined below. Debye = 9N Av k B ( knowledge of the specific heat and the sound velocity allows to infer the mean free path from the thermal conductivity.APPENDIX E: FREQUENCY DEPENDENT SPECULARITY PARAMETER Instead of using a constant specularity parameter as explained in the main text, the introduced formalism can be extended to account for its frequency dependence by using the same roughness distribution (ηBy relating the wavelength to the frequency as λ = 2π ω , one can therefore introduce a frequencydependent specularity parameter (ω) in the effective boundary length   eff (ω) = 1+(ω) 1−(ω)  ′ entering in Eq. (1) which is then integrated numerically.APPENDIX F: CALCULATION OF THE CASIMIR LIMIT FROM DOMAINDISTRIBUTIONSThe Casimir regime assumes that incident phonons which are absorbed by the various surfaces of the sample are then re-emitted with the equilibrium distribution corresponding to the local temperature.Since the emitted energy varies as T4 , the temperature can be expanded as T4 (z) = T 4 (0) + 4T 3 (0) z (dT/dz) by assuming a constant temperature gradient in the z direction, with z=0 the center of the sample.By dividing the emitted energy flux by the cross-section S and (dT/dz), one can then write the thermal conductivity   in the Casimir limit regime.
Casimir regime is recovered by relating the latter to the configuration factors Fzi.

3 2π(
+   +   ) by introducing the configuration factors associated to the three kinds of surfaces.Fzz is in particular defined below with the width a (along x), the thickness b (y), the z-length c, with the notation =2a/c and =2b/c.

FIG. 1 . 3 𝐶
FIG. 1. Optical images in reflection of selected single crystals of LaAlO3.(a) after annealing, (b,c) as-received, (e,f) after quenching.Panels are oriented such that the heat flow in thermal conductivity measurements goes from the left to the right, as indicated by black arrows.(d) Schematics of the domain structure showing two domains, where ellipses indicate the orientation of their axes of compression (minor axis) and extension (major axis), projected into the plane of the schematic.Scale bars correspond to 100 µm.(g) Histogram of the number of domains with given sizes, measured from optical images.

FIG. 3 .
FIG. 3. (a,b) Thermal conductance as a function of temperature and (c,d) mode transmittance as a function of phonon frequency in the harmonic approximation, as simulated using non-equilibrium Green's function techniques.The black lines correspond to LaAlO3 in the monodomain state, and blue (respectively red) lines in (a) and (c) [respectively (b) and (d)] correspond to LaAlO3 with domain walls orthogonal (respectively parallel) to the heat flow.

TABLE 2 .
Boundary lengths.Average domain size measured on optical images.Boundary lengths with constant specularity parameter p for LB (0=7 nm) or with a frequency dependent one p() for LB' (0=4 nm).Casimir limit LC inferred from the optically determined domain distribution and Fig.4(d).