Carmiña Gartner

Betty Lucy López

Ligia Sierra

Robert Graf

Hans W Spiess

Marianne Gaborieau
email: m.gaborieau@uws.edu.au

Interplay between structure and dynamics in chitosan films investigated with solid-state NMR, dynamic mechanical analysis and X-ray diffraction

Keywords: chitosan, solid-state NMR, local dynamics, interactions, dynamic mechanical analysis, X-ray diffraction

Modern solid-state NMR techniques, combined with X-ray diffraction, revealed the molecular origin of the difference in mechanical properties of self-associated chitosan films. Films cast from acidic aqueous solutions were compared before and after neutralization, and the role of the counterion (acetate vs Cl -) was investigated. There is a competition between local structure and long-range order. Hydrogen bonding gives good mechanical strength to neutralized films, which lack long-range organization. The long-range structure is better defined in films cast from acidic solutions, in which strong electrostatic interactions cause rotational distortion around the chitosan chains. Plasticization by acetate counterions enhances long-range molecular organization and film flexibility. In contrast, Cl -counterions act as a defect and impairs the long range organization by immobilizing hydration water. Molecular motion and proton exchange are restricted, resulting in brittle films despite the high moisture content.

INTRODUCTION

Chitin, the second most abundant natural polysaccharide, is synthesized by a number of living organisms and is the major structural component in the exoskeleton of arthropods and cell walls of fungi and yeast. The main commercial sources of chitin are crab and shrimp shells, which are abundantly supplied as waste products of the seafood industry 1 . Because chitin is not readily dissolved in common solvents, it is often converted to its partially deacetylated derivative, chitosan, a copolymer of glucosamine and N-acetylglucosamine (Figure 1). Chitin and chitosan are biocompatible, biodegradable, and non-toxic, and have anti-microbial and hydrating properties. Chitosan is easily processed into films, gels, nanofibers, beads, microparticles, nanoparticles and scaffolds and thus widely used for biomedical and pharmaceutical applications, tissue engineering, controlled drug delivery, biotechnology, and in the food industry 2-4 . Chitosan can be formed into hydrogels or films using an external crosslinker, but it can also undergo a self-crosslinking process. The organization of the chitosan self crosslinked structure is due to inter-and intramolecular associations via hydrogen bonding or hydrophobic interactions and has been the subject of numerous studies [5][6][7][8][9][10][11] . Nevertheless the self-crosslinking mechanism remains incompletely solved.

The intra-and intermolecular interactions also govern the macroscopic thermal and mechanical properties of chitosan films, therefore it is important to elucidate them. Due to the hygroscopic nature of chitosan 12 , it always contains residual water (more strongly bound to hydroxyl than to amine groups 13 ), and its characterization must be performed at defined moisture content [14][15] . The residual water impairs the determination of the glass-transition and secondary relaxation phenomena in this semi-crystalline polysaccharide, explaining the scatter of published data obtained by differential scanning calorimetry (DSC), dynamic mechanical analysis (DMA) and dielectric relaxation spectroscopy (DRS). [16][17][18][19][20] This variation is related to the complex structure in chitosan via the local dynamics, which is still to be elucidated.

NMR is one of the most powerful techniques to access to site-specific structure and dynamics of complex molecular architectures. 1 H NMR, 13 C NMR, and 15 N NMR in solution and in the solid have being widely used to quantify the degree and sequences of chitosan acetylation [21][22][23][24][25] . Other less common studies have also been reported such as: interactions of chitosan with crosslinking agents 26 , water 27 and plasticizers 28 , studies of new chitosan copolymers structures 29 , and conformational studies 30 . However, the local dynamics and interactions in self-crosslinked chitosan films have not yet been reported by solid-state NMR. Solid-state NMR experiments performed under fast magic angle spinning (MAS) provide new and versatile spectroscopic tools, taking advantage of the significant simplification of the dipolar-coupled multispin network due to its decoupling performance. Recoupling techniques that involve appropriate pulse sequences have been applied along with fast MAS to obtain structural and dynamics information inherent to 1 H-1 H and 1 H-13 C dipolar couplings, recovering relevant interactions selectively [START_REF] Schnell | 1H Solid-State NMR of Supramolecular Systems[END_REF] . Specifically, the rotational-echo double-resonance (REDOR) is a recoupling technique routinely used to study biological systems, e.g. to obtain internuclear distance constraints between labeled spin pairs in biomolecules. [START_REF] Gullion | [END_REF] Hydrogen bonding structures in supramolecular systems have been studied through the recoupling of homonuclear dipole-dipole interactions, using the so-called back-toback pulse sequence for 1 H double quantum (DQ) spectroscopy experiments 33 . In such a twodimensional experiment, double quantum coherences due to pairs of dipolar coupled protons are correlated with single quantum coherences, resulting in correlation peaks characteristic for 1 H-1 H pairs. A variant of the REDOR sequence can be implemented to recouple coherent polarization transfer for heteronuclear single-quantum correlation experiment to identify the directly bonded carbon resonances 34 .

Even though the overall mechanism occurring in the formation of self-assembling chitosan films and gels is partially identified, little is known about the interactions that remain after their formation, and about their influence on the physicochemical properties of the films. We aim at gaining insight into the complex self-association network in chitosan films cast from dilute solutions of a weak and a strong acid at approximately the same concentration. Non-neutralized and neutralized films exhibit very different mechanical properties. Their systematic study with modern NMR techniques will probe their local structure and dynamics. Complemented by the mechanical properties investigated by DMA and the long-range organization studied by XRD, it will shed light on the molecular interactions in solid chitosan films and the role of the possible counterion.

EXPERIMENTAL SECTION

2.1. Sample Preparation. Chitosan was purchased from Aldrich (medium molecular weight, catalog number 448877, batch 09303PE). Its degree of acetylation DA (see Figure 1) of 31% and its molecular weight of 240 kDa were previously measured by 13 C solid-state NMR and viscosimetry respectively 35 .

It was used to prepare 3 w/v % chitosan solutions (above the entanglement concentration, experimentally found as 1.2 w/v % 36 ) using two types of dilute acids: 0.17 M acetic acid (AcOH) and 0.12 M hydrochloric acid (HCl). Films were cast from these solutions and dried at room temperature.

The neutralization of the films was performed by immersing the dried films in a 1 w/v % solution of NaOH, and then thoroughly washing with distilled water. The dried films were equilibrated under 44% controlled relative humidity in a desiccator with a saturated solution of K 2 CO 3 . Films are denoted: F HCl (film cast from HCl solution), F AcOH (film cast from acetic acid solution), FN HCl (film cast from HCl solution then neutralized) and FN AcOH (film cast from acetic acid solution then neutralized).

Note that films were also cast from more dilute solutions of chitosan (1 w/v%, below the entanglement concentration, with all other parameters unchanged); the NMR spectra recorded show no significant difference with the films prepared from 3 w/v% chitosan solution, and will only be shown in supporting information.

Thermogravimetric Analysis (TGA).

TGA measurements were carried out on a TGA 500 (TA Instruments). The moisture content was measured by heating the samples to 120 °C under nitrogen atmosphere using a ramp of 2 °C•min -1 ; at this temperature an isotherm was kept for 60 min. The percentage of moisture was calculated from the weight lost at the end of the isotherm. To study the films degradation with temperature, samples were run in air from 50 to 800 °C, using a ramp of 10 °C•min -1 . 15 N cross polarization under magic-angle spinning (CP-MAS) NMR experiments were performed on a Bruker Avance II console operating at 300 MHz 1 H Larmor frequency with a commercial double resonance probe supporting zirconia MAS rotors with 4 mm outer diameter at a MAS spinning frequency of 10 kHz and 62.5 RF nutation frequency (5 µs 90 º pulse length), 5 ms contact time and 5 s relaxation delay. 13 C CP-MAS NMR experiments were performed on a Bruker Avance console operating at 500 MHz 1 H Larmor frequency with a commercial double resonance probe supporting zirconia MAS rotors with 2.5 mm outer diameter at a MAS spinning frequency of 25 kHz and 125 RF nutation frequency (2.5 µs 90 º pulse length), 1 ms contact time and 5 s relaxation delay. 1 H MAS NMR experiments were performed on a Bruker Avance III console operating at 850 MHz 1 H Larmor frequency with a commercial double resonance probe supporting zirconia MAS rotors with 2.5 mm outer diameter at a MAS spinning frequency of 30 kHz and 125 RF nutation frequency (2.5 µs 90 º pulse length). 1 H double quantum correlation spectra as well as REDOR-based heteronuclear ( 1 H-13 C) correlation spectra were recorded at 50 kHz MAS and 700 MHz 1 H Larmor frequency on a Bruker Avance spectrometer using a Bruker 1.3 mm double resonance MAS probe with 140 kHz RF nutation frequency (1.8 µs 90 º pulse length). Double quantum excitation was achieved using the Back-to-Back recoupling pulse sequence with 40 s recoupling time, whereas the REPT-HSQC sequence has been applied for 1 H-13 C heteronuclear correlation experiments with 80 s recoupling time.

Nuclear Magnetic Resonance (NMR).

Dynamic Mechanical Analysis (DMA).

The mechanical properties were measured with a TA instruments DMA Q800, using the tension clamp for rectangular film and a temperature sweeping from 40 to 200°C, under the following conditions: heating ramp 3 °C•min -1 , initial strain 0.01%, preload force 0.3 N and frequency 1 Hz. Three measurements were done for each film.

X-Ray Diffraction (XRD).

Diffractograms were obtained in a Phillips PANalytical X'PertPro MPD, with a wavelength of 1.542Å (Cu K 1 1 radiation), a 2 range from 5 to 50° and a resolution of 0.013°. Diffractograms were normalized with sample amount.

RESULTS AND DISCUSSION

3.1. Moisture Content. The moisture content was evaluated by the decrease of sample weight at the end of an isotherm at 120 °C. The stabilization of weight lost within the first 60 minutes of the isotherm confirmed that this program was accurate in terms of moisture determination. 120 °C was low enough to prevent thermal degradation of the materials. The moisture content was 16 % for F AcOH and 18 % for F HCl , while for the neutralized films was lower: 13 % for the FN AcOH and 12 % for FN HCl . The moisture content of the chitosan powder was 10 %. Although the moisture content is similar, it is slightly higher for the protonated films, which may be explained by the hydrophilicity of the NH 3 + groups. The quantified moisture contents are considered to be low moisture content 14 .

Solid-state NMR Results.

3.2.1. 15 N and 13 C NMR. In diluted acidic media, chitosan contains glucosamine units both in the protonated and neutral amine forms. With a weak acid, such as acetic acid, a fraction of the amine groups will be protonated, whereas all amine groups will be protonated with a strong acid, such as hydrochloric acid. The films cast from chitosan solutions could thus contain charged and neutral amine groups. 15 N CP-MAS spectra allow for a direct detection of each monomer unit: the acetylated amine at about 100 ppm, the amine between 0 and 15 ppm (Figure 2a). It also discriminates between charged and neutral amine groups, detected respectively at 15 and 0 ppm. It is experimentally observed that the films directly cast from chitosan solutions in HCl are charged, while the original powder and the films suspended in dilute NaOH solutions and rinsed are neutral.

13 C CP-MAS NMR spectra show 'snapshots' of the distribution of chemical shifts present in the sample. 13 C CP-MAS spectra were recorded on the original powder, films prepared with acetic acid and HCl, and neutralized films (see Figure 2b for representative spectra and supporting information for all spectra). Powder and neutral films all exhibit the same 13 C CP-MAS spectrum. A full chemical shift assignment can be found in the literature. 21

One-dimensional 1 H NMR.

1 H MAS spectra were recorded at different temperatures on the original powder, films prepared with acetic acid and HCl and neutralized films at high magnetic field under fast MAS (see Figure 3 for representative spectra and supporting information for all spectra). At room temperature (Figure 3a), all samples apart from F HCl exhibit broad signals due to a large range of local packing arrangements. In FN AcOH , the acetate was washed away. In F HCl , the much more narrow linewidths clearly reveal an NH 3 + signal at 8.5 ppm (and a methyl signal at 2.2 ppm).

The temperature behavior of the powder and the neutralized films are the same (see Figure 3b and supporting information). The only significant change observed upon heating is a reversible broadening of the signal originally at 4.3 ppm, accompanied by its reversible shift to 3.9 ppm. This signal originates in the residual water and the OH groups of chitosan (as confirmed by its disappearance in a chitosan powder deuterated by suspension in D 2 O, see supporting information). Note that this signal is weaker after cooling down, probably due to partial evaporation of the residual water. The rather unusual line broadening upon heating indicates that the water molecules and the OH groups in the sample are in the slow exchange limit with another site in the sample, most likely the NH 2 groups. Here the slow exchange limit means: there is basically no 1 H exchange taking place at room temperature on the time scale of the experiment, and local molecular motions as well as chemical exchange of protons start on that time scale upon heating.

The narrow signals of sample F HCl (Figure 3c) could originate either in a strong local mobility (resulting in limited line broadening), or in a well-organized local structure (resulting in a limited distribution of chemical shifts). As above, the broadening of this signal upon heating indicates limited mobility at room temperature, and the onset of local motion or chemical exchange at higher temperatures. The higher resolution obtained for this sample allows observing that only two signals broaden upon heating, to a similar extent: the signal of residual water and OH groups at 4.7 ppm on one hand, and the NH 3 + signal at 8.3 ppm on the other hand. Furthermore, the NH 3 + signal shifts upfield towards the signal of residual water and OH groups (see insert), while the signal of residual water and OH groups shifts upfield (as a result of a shift downfield towards the NH 3 + signal and a stronger shift upfield due to the temperature increase). This supports the interpretation of a chemical exchange taking place between these. The narrow NH 3 + signal observed for F HCl compared to NH 2 in other samples (around 8 ppm) is due to a greater distribution of chemical shifts for NH 2 due to a large variety of hydrogen bonding structures in other samples, as opposed to more similar electrostatic interactions through NH 3 + in F HCl . Interestingly, F HCl is the film with the highest moisture content but not the film with the highest mobility, although water could be expected to plasticize the film. This indicates that a significant part of the residual water is immobilized in this particular sample. As it is not the case to that extent in the other samples, it may be bound to the Cl -counterion.

In F AcOH , the signals of the acetate (5.3 ppm and 2.1 ppm) become narrower upon heating (Figure 3d), as is expected when the local molecular motion becomes faster. Thus there is a significantly stronger mobility in F AcOH than in all the other samples. It can be concluded that more than the charge, the counterion plays a significant role in the local dynamics in chitosan films. The acetate counterion's charge is delocalized over the carboxylic acid group, the ion is larger and has more degrees of freedom (thus can act as a plasticizer).

3.2.4.

Two dimensional 1 H-13 C NMR. The chitosan powder, F AcOH and FN HCl show the strongest dipolar interaction between 75 and 80 ppm, which corresponds to C5, and a weaker dipolar interaction between 50 and 60 ppm, corresponding to C2 carbons (Figure 5). The relative rigidity of these two carbons is mostly intrinsic of the chitosan structure, as they indicate the rigidity of the glucose ring.

Furthermore, earlier X-ray diffractions studies [39][40] indicated strong intramolecular hydrogen bonding in chitosan through oxygen bound to C5 and the amino group bound to C2 carbons. In F HCl , there is an additional strong signal around 70-75 ppm corresponding to C3, and an increase in the interaction for the C2 carbons (near 55 ppm). This indicates the rigidity around the NH 3 + group and the hydration shell of its Cl -counterion.

Figure 5. 1 H-13 C HSQC REPT correlation spectra for the chitosan films and powder, recorded with a recoupling time of 80 s.

DMA Results

. F HCl exhibits different properties from the other films, as it is visually more brittle.

All films were investigated with DMA at different temperatures to investigate their viscoelastic behavior in more detail. Storage moduli E' are summarized in Table 1. E' is higher for F HCl than for F AcOH , while FN HCl and FN AcOH exhibit an intermediate behavior. This confirms the visual stiffness of F HCl compared to the other films, due to the strong electrostatic interactions around the NH 3 + groups and their Cl -counterions, immobilizing hydration water. In FN AcOH and FN HCl the chain-chain interactions lead to an intermediate stiffness, while in F AcOH the acetate counterion acts as a plasticizer, imparting more elasticity and flexibility to this film. and FN HCl can be explained by its higher moisture content, and the rapid chemical exchange between NH 2 and NH 3 + groups. The film with the highest moisture content, F HCl , however, does not follow this trend. We have demonstrated above that the water molecules in F HCl are tightly bound to the chitosan and its Cl -counterion, thus the loss of water by this film is very limited around 100 °C. Furthermore, all amine groups are protonated in F HCl , strongly limiting the possibility of chemical exchange between NH 2 and NH 3 + groups. The strong attenuation of the first relaxation for F HCl thus indicates that water and H + are more strongly bound in this film.

The transition around 170 °C is only seen for the charged films. It is pronounced for F AcOH , and it is only beginning for F HCl when the film ruptures. To the best of our knowledge, this transition was not reported before. It could be linked to hydrogen bonding (which exist at high temperatures, provided the samples are rigid enough) or to electrostatic interactions due to the presence of the NH 3 + groups. Work is in progress to record NMR data in this temperature range to gain molecular information. The rupture of F HCl at a temperature at which the other films are stable correlates with the onset of its decomposition at ca 400°C (compared to ca 500°C for the other films) as observed by TGA.

These differences in DMA behavior lead us to conclude that the NH 3 + and NH 2 side groups have a drastic influence on the local dynamics. They determine the relative amount of the more strongly bound water and the nature of the hydrogen bonding network formed between chitosan moieties and water. There is an apparent contradiction in stating that charged films tend towards a much more regular long-range order, but that the least charged one is the best organized. It can be rationalized by taking into account that the strong rigidity of F HCl may prevent the structure to reorganize, which is needed to accomplish long-range order.

Note that the variations in 13 C CP-MAS spectra with the acid solution from which the films were cast (HCl and AcOH) had been previously assigned to different helical structures. 30 However, they exhibit the same X-ray diffractogram, thus the same structure. We have shown here that this is due to restricted mobility in F HCl , hindering the long-range organization of the film in the structure it tends to. This is in agreement with the previously observed increase in the intensity of the X-ray reflections upon annealing 30 , which we can now assign to a better long-range organization owing to increase local mobility.

CONCLUSIONS

Modern solid-state NMR techniques, combined with X-ray diffraction, revealed the molecular origin of the difference observed in macroscopic properties of chitosan films. Films cast from acidic aqueous solutions were compared before and after neutralization, and the role of the counterion (acetate vs Cl -) was investigated. There is a competition between local structure and long-range order in self-associated chitosan films. Neutralized films exhibit a molecular structure similar to that of the original chitosan powder, not organized on the long range, mediated through intra and intermolecular hydrogen bonding.

These films show good mechanical strength and elasticity. The chitosan is fully charged in films cast from HCl, while it is only partially charged in films cast from acetic acid. The long-range structure is better defined in films cast from acid solutions, in which the presence of charges induces strong electrostatic interactions and rotational distortion around the main chain, stronger in the case of HCl.

The large acetate counterion with a delocalized charge plasticizes the films enabling a good long-range molecular organization, as well as increased film flexibility. The Cl -counterion acts as a defect and impairs the long range organization by strongly immobilizing water in its hydration shell. Both the molecular motion and the proton exchange are very restricted, resulting in the film being stiff and brittle, although it exhibits the highest moisture content of all films.

The changes in local conformation and dynamics introduced by the charge in the chitosan chain are reversible. The increase in mobility and in proton exchange in the films upon heating is also reversible, even after partial moisture evaporation. The chitosan concentration in the initial solution does not affect the structure and dynamics in the films, as shown by the similar behavior of films cast from chitosan solution below and above the entanglement concentration. A relaxation process was detected at about 170 ºC for charged films, resulting in film rupture in the case of HCl. It could be linked to hydrogen bonding or to electrostatic interactions; NMR data are needed in this temperature range to reach any further conclusions.

The innovative methodology developed in this work can be extended to the investigation of the effect of other parameters on the interactions in self-association chitosan systems. Films obtained at higher degrees of acetylation, hydrogels obtained at higher moisture content, should exhibit different local structure and dynamics, which are expected to change the macroscopic properties. It will be the subject of further studies.

Figure 1 .

 1 Figure 1. Chitosan molecular structure. DA is the degree of acetylation.

  Two more signals are present in films prepared with acetic acid (C9 at 25 ppm, C8 at 180 ppm) due to residual acetic acid. Charged films exhibit some signals shift as compared to neutral ones, and these shifts are stronger for F HCl than for F AcOH . The signal of C1 originally at 105 ppm shifts to 100 ppm then 98 ppm, the signal of C4 originally at 83 ppm shifts to 82 then 81 ppm, the signal of C3 originally at ca 74 ppm shifts to 72.5 ppm then 71.5 ppm. Such shifts were already observed30 and assigned to different helical structures in the solid. It will be discussed below (section 3.4.) that both F HCl and F AcOH exhibit the same XRD diffractogram (albeit with different intensities), thus it can hardly be concluded the signals shifts originate in such a radical structural change. To check if the signals shifts were not merely a result of the various degrees of protonation of the chitosan, chemical shifts observed in 13 C CP-MAS were compared to published chemical shift values measured in solution. Chitosan is not soluble in water unless it is charged, thus only spectra recorded in acidic solutions can be compared. Signals of C1, C4, C3 and C2 were observed at 98, 78, 71 and 56 respectively in D 2 O with AcOD 37 , while they were observed at 100, 80, 73 and 58 in D 2 O with DCl38 . The signal of C2 is not expected to shift with pH (it does not shift in 13 C CP-MAS), thus the published spectra are simply calibrated differently, and the signals of C1, C4 and C3 in solution do not depend on whether the acid is strong (DCl) or weak (AcOD). Therefore the shifts observed with 13 C CP-MAS are due to packing effects. The strongest effect is observed on the signals of C1 and C4, which are the carbons linking adjacent glucosamine rings with a glucosidic bond. Thus it can be explained by a rotational distortion around the main chain, the presence of the charge introducing strong electrostatic interactions which could hinder the packing of the chain following the usual hydrogen bonding pattern.The more charged film, F HCl , thus exhibits stronger rotational distortions and stronger signals shifts than F AcOH . The local mobility in films will be examined below. A smaller shift is observed for C3, which is spatially close to the charge.15 N and 13 C NMR measurements indicated charge and conformational changes in the films, and opened questions on the local molecular dynamics. 1 H NMR was employed to investigate the dynamics in the films and the different association mechanisms.

Figure 2 .

 2 Figure 2. Solid-state NMR spectra of chitosan powder and films: 15 N CP-MAS (a) and 13 C CP-MAS (b). Vertical dashed lines in (b) are a guide to the eye.

Figure 3 .

 3 Figure 3. 1 H MAS NMR spectra of various samples at 300 K(a); temperature-dependent 1 H MAS NMR spectra recorded for the powder (b), F HCl (c) and F AcOH (d). For temperature-dependent spectra, the spectra were acquired in the order of the legend (top to bottom). The insert in (c) shows a zoom in the NH 3 + signals recorded upon heating.

Figure 4 .

 4 Figure 4. Two-dimensional 1 H DQ MAS NMR correlation spectra for the chitosan films and powder recorded with the back-to-back sequence and a recoupling time of 40 s.

Figure 6 .

 6 Figure 6. Variation of the loss tangent with temperature for the chitosan films.

Figure 7 .

 7 Figure 7. X-ray diffractograms of chitosan films.

Table 1 .

 1 Storage moduli E' measured at different temperatures by DMA (average of threeThe temperature dependence is shown for the loss tangents (tan ) on Figure6. The loss factor for F AcOH has two maxima at ~ 95 °C and ~ 172 °C, which can be related to conformational relaxations. F HCl only exhibits a very weak transition at ~ 100 °C and ruptures at ~ 170 °C, while for the neutralized films FN HCl and FN AcOH there are two overlapping broad relaxations with maxima at ~ 97 °C and 104 °C. The relaxation observed at ~ 97 ºC has been assigned in previous dielectric spectroscopy studies to water desorption 12 and to conduction via H + migration20 . Its stronger intensity in F AcOH than in FN AcOH

	measurements).
	T (°C)			E' (GPa)
		F HCl	FN HCl	F AcOH	FN AcOH
	50	4.1 ± 0.2	2.6 ± 0.6	1.7 ± 0.2	3.4 ± 0.4
	100	4.2 ± 0.4	2.3 ± 0.5	1.8 ± 0.2	2.9 ± 0.3
	195	-	2.2 ± 0.4	2.1 ± 0.2	2.5 ± 0.2

H-13 C HSQC REPT correlation spectrum of a chitosan film prepared from 1% chitosan in HCl solution.

ACKNOWLEDGMENTS.

We are grateful to the Deutscher Akademischer Austausch Dienst (DAAD), the University of Antioquia and the Max Planck Society for financial support.

Two-dimensional solid-state NMR spectra

Two-dimensional 1 H DQ MAS NMR correlation spectrum of a chitosan film prepared from 1% chitosan in HCl solution, recorded with the back-to-back sequence and a recoupling time of 40 s.

SUPPORTING INFORMATION AVAILABLE: solid-state NMR spectra, TGA and DMA thermograms, X-ray diffractogram. This material is available free of charge via the Internet at http://pubs.acs.org.

(1) Rinaudo, M. Prog. Polym. Sci. 2006, 31, 603-632.

(2) Jayakumar, R.; Prabaharan, M.; Nair, S. V.; Tamura, H. Biotech. Adv. 2010, 28, 142-150.

Supporting information for

Interplay between structure and dynamics in chitosan films investigated with solid-state NMR, dynamic mechanical analysis and X-ray diffraction Carmiña Gartner, Betty Lucy López, Ligia Sierra, Robert Graf, Hans W. Spiess, and Marianne Gaborieau TGA thermograms 15

N CP-MAS NMR spectra

The chitosan concentration (1% or 3%) does not have an influence on the 15 N CP-MAS NMR spectrum.

S-2 13 C CP-MAS NMR spectra

Powder and neutral films all show the same spectrum Films prepared in acetic acid show 2 more signals (CH 3 at 25 ppm, C=0 at 180 ppm) due to residual acetic acid.

Films which are charged show some signals shift as compared to neutral ones, and these shifts are stronger for films prepared in HCl than in AcOH.

HMAS NMR spectra

All samples measured at 300, 340, 380, 340, 300 K in that order (blue, red, green, purple then yellow).

Temperature-dependent 1 H MAS spectra of chitosan powder. Temperature-dependent 1 H MAS spectra of a film prepared from 1% deacetylated chitosan in HCl solution.

Influence of the pH on the mechanical behavior of the chitosan cast films

Two films were cast from solutions of 3% w/v chitosan in HCl and acetic acid at the same pH value of 4.2. The mechanical performance of these dry films were measured by DMA (using the same conditions that were reported), to verify if the behavior previously observed depended on the pH of the solution. It is observed in the figure below that the mechanical behavior of both films does not change from what was reported before and that it is independent from the pH. There are two relaxations for the film from AcOH and nearly no observed relaxations for the film from HCl.

Films cast from HCl solutions at pH of 2.0 and 1.5, were also measured by DMA (films were not formed at lower pHs). These films displayed the same mechanical behavior trend. The only difference was a slight change on the onset of the degradation temperature, which is higher at higher pH. These values were 152 °C for pH 1.5, 166 °C for pH 2.0 and 182 °C for pH 4.2.

It can be concluded that of dissociation of the acid as well as the nature of the counterion play a major role in the mechanical properties of chitosan cast films.

Comparison of the DMA thermograms for the films cast from HCl and AcOH solutions at pH 4.2.

Comparison of the DMA thermograms for the films cast from HCl solutions at different pH.


domain: spi math

text:
Grigor Nika
email: grigor.nika@kau.seadrian.muntean@kau.se

Adrian Muntean

Hypertemperature effects in heterogeneous media and thermal flux at small-length scales

Keywords: MSC 2020 74Q05, 35B27, 35G05, 35G15, 35J58, 35Q79 correctors, scale-size thermal effects, generalized Fourier's law, microstructure

We propose an enriched microscopic heat conduction model that can account for size effects in heterogeneous media. Benefiting of physically relevant scaling arguments, we improve the regularity of the corrector in the classical problem of periodic homogenization of linear elliptic equations in the three-dimensional setting and, while doing so, we clarify the intimate role correctors play in measuring the difference between the heterogeneous solution (microscopic) and the homogenized solution (macroscopic). Moreover, if the data are of the form f = div F F F with F F F ∈ L 3 (Ω, R 3 ), then we recover the classical corrector convergence theorem.

Introduction

Analysis of correctors in homogenization theory for second-order elliptic equations with highly oscillatory coefficients provides an important link between the microscopic and macroscopic problem since it measures the difference between the heterogeneous solution and the homogenized solution. Furthermore, it provides key information needed for proving the convergence of multiscale approximation schemes based on, e.g., heterogeneous multiscale method (HMM) or multiscale FEM (MsFEM). Simultaneously, from a modeling perspective, correctors point out some of the limitations of such second-order elliptic A c c e p t e d V e r s i o n systems due to their lack in accounting for scale-size effects in e.g. micro-heterogeneous bodies. The prototypical and landmark example for the homogenization of second-order elliptic equations is Fourier's law of heat conduction. In its simplest form, Fourier's law relates the heat flux q q q as a linear function of the temperature gradient, that is:

q q q = -κ∇u, (1.1)
where u is the absolute temperature and κ > 0 the thermal conductivity that depends on properties of the material. In general, the proportionality coefficient κ (either scalar or tensor) may depend on temperature, space and/or time variables, or on other parameters, but often varies so little in cases of interest that it is reasonable to neglect this variation.

The theory of periodic homogenization led to a deeper exploration of Fourier's heat conduction law as it pertains to heterogeneous periodic material with different conductivities. Effective heat fluxes were derived taking into account microstructure morphology and volume fraction. The classical problem in the periodic homogenization for stationary heat conduction states: Find u ε ∈ H 1 0 (Ω) satisfying,

-div K( x x x ε ).∇u ε = f in Ω, u ε = 0 on ∂Ω, (1.2)
where K(y y y) ∈ L ∞ (Y, R 3×3 ) is uniformly elliptic, symmetric, and Y -periodic with Y = [0, 1) 3 , and explore what happens with the situation when ε → 0. If Ω ⊂ R 3 is uniformly Lipschitz open set, then there exists a unique solution u ε to (1.2) such that it converges weakly to a function u in H 1 0 (Ω), where the function u ∈ H 1 0 (Ω) is the unique solution to,

-div K eff .∇u = f in Ω, u = 0 on ∂Ω, (1.3)
with K eff ℓj := 3 i=1 Y K ℓi (y y y)(δ ji -∂ y i w j ) dy y y and w j ∈ {H 1 per (Y ) | Y w j dy y y = 0} solution to, -div y (K(y y y).(∇ y w j -e e e j )) = 0 in Y, w j is Y -periodic.

(1.4)

The convergence of u ε to u in H 1 0 (Ω) is only weak. If one seeks to improve the convergence then, usually, a corrector type term is introduced,

u ε -u -εû(•, • ε ) → 0 in H 1 (Ω), (1.5)
where û(x x x, y y y):= -N k=1 w k (y y y) ∂u ∂x k (x x x) is the correction term. In order for the expression to belong in H 1 (Ω) it is required that w k ∈ W 1,∞ (Y ) and u ∈ H 2 (Ω) (see [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF]Sect. 5, A c c e p t e d V e r s i o n pg. 33]). With the introduction of two-scale convergence in [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF] (see also [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] and [START_REF] Lukkassen | Two-scale convergence[END_REF]) a rigorous justification was provided for the multiple-scale method and the corrector result (in N -dimensional space) in (1.5) was made rigorous through the following theorem: Proposition 1.1. [2, Thm 2.6], [10, Prop. 9.12], [START_REF] Lukkassen | Two-scale convergence[END_REF]Thm. 22] Let ũ1 ∈ L s (Ω) and û be given by,

û(x x x, y y y) = - N k=1 w k (y y y) ∂u ∂x k (x x x) + ũ1 (x x x), (1.6)
and suppose that ∇ y w k ∈ L r (Y, R N ), k = 1, . . . , N and ∇u ∈ L s (Ω, R N ) with 1 < r, s < ∞ and such that,

1 r + 1 s = 1 2 . (1.7)
Then,

∇u ε -∇u -∇ y û(•, • ε ) → 0 in L 2 (Ω, R N ). (1.8)
As one can immediately observe, the main unease with the above proposition is the higher integrability required on both local and homogenized solutions. However, with the introduction of periodic unfolding operators [START_REF] Ciorȃnescu | Éclatement périodique et homogénéisation[END_REF][START_REF] Ciorȃnescu | The periodic unfolding method in homogenization[END_REF], one can obtain a more general corrector result without requiring any regularity assumption on the cell function w k , stating,

∇u ε -∇u - N k=1 Q ε ( ∂u ∂x k )∇ y w k ( • ε ) → 0 in L 2 (Ω, R N ), (1.9)
where Q ε is the scale-splitting operator defined in [START_REF] Griso | Error estimates and unfolding for periodic homogenization[END_REF]; see the elegant proof in [START_REF] Ciorȃnescu | The Periodic Unfolding Method. Theory and Applications to Partial Differential Problems[END_REF]. Moreover, based on the scale-splitting operator, upper bound estimates on the convergence rate in terms of ε were obtained in [START_REF] Griso | Error estimates and unfolding for periodic homogenization[END_REF] (with some additional regularity assumptions on the homogenized solution). Furthermore, the upper bound estimates on the convergence rate can be made tighter by using boundary layer correctors (see, e.g., [START_REF] Moskow | First-order corrections to the homogenised eigenvalues of a periodic composite medium. a convergence proof[END_REF], [START_REF] Allaire | Boundary layer tails in periodic homogenization[END_REF], [START_REF] Onofrei | Error estimates for periodic homogenization with non-smooth coefficients[END_REF]). Therefore, it seems that in order to lift the (restrictive) regularity conditions in Proposition 1.1 one must have knowledge of the operator Q ε .

The widely accepted understanding of Fourier's law of heat conduction is that it represents a limiting case approximation of a more comprehensive, potentially nonlinear, constitutive law for the heat flux, which may depend on higher-order gradients (see, e.g., [START_REF] Coleman | Thermodynamics and departures from Fourier's law of heat conduction[END_REF], [START_REF] Christov | On a higher-gradient generalization of Fourier's law of heat conduction[END_REF], [START_REF] Pan | A generalized heat conduction model in rarefied gas[END_REF]). For instance, the temperature of a rarefied gas at the slip regime, namely when 0.001 < Kn < 0.3, where Kn is the non-dimensional Knudsen number, deviates from Fourier's law of heat conduction (see, e.g., [START_REF] Pan | A generalized heat conduction model in rarefied gas[END_REF]). Moreover, in the same article a generalized heat conduction model, from a phenomenological point of view, was postulated under the assumption that the gas is isotropic. The authors' speculation was that the heat flux in a rarefied gas in the slip regime, depends linearly on the temperature gradient but also on higher-order temperature derivatives,

A c c e p t e d V e r s i o n q q q = K.∇u + L:∇∇u + M . . . ∇∇∇u, (1.10) where K is the classical second order heat conduction tensor1 , while L and M are third and fourth order material parameter tensors, respectively (the notation ., :, and . . . denote first, second, and third order contractions). The structure of such material parameter tensors can be determined by considering proper orthogonal transformations (compare [START_REF] Suiker | Application of higher-order tensor theory for formulating enhanced continuum models[END_REF]). Moreover, if we request that the material be centro-symmetric (meaning it is invariant under reflections), then the tensor L ≡ 0 (cf. [START_REF] Suiker | Application of higher-order tensor theory for formulating enhanced continuum models[END_REF]). Furthermore, if the material is isotropic then reflections are not relevant since the only isotropic third order tensor is the Levi-Civita tensor which is anti-symmetric and, hence, L:∇∇u ≡ 0 since it would be the product of a symmetric and an antisymmetric matrix which is zero. Additionally, the tensor M has a physical meaning -it is often referred to as the spatial retardation (see the context of [START_REF] Christov | On a higher-gradient generalization of Fourier's law of heat conduction[END_REF]).

Henceforth, we assume isotropy of the medium. Thus, the higher-gradient heat flux can reduce further,

q q q = K.∇u + M . . . ∇∇∇u, (1.11)
where K ij :=ηδ ij and M ijkl :=ℓ 1 δ ik δ jl + ℓ 2 δ il δ jk + ℓ 3 δ ij δ kl for scalars η and ℓ i , i = 1, 2, 3 that will be assumed to be constant or piece-wise constant with the heat conduction inequality asserting, q q q.∇u ≤ 0, (1.12) for all temperature fields.

Non-classical laws of Fourier's heat conduction have, for many years now, attracted considerable attention from the theoretical mechanics community (see, e.g., [START_REF] Aifantis | Internal length gradient (ILG) material mechanics across scales and disciplines[END_REF], [START_REF] Berezovski | Internal variables representation of generalized heat equations[END_REF], [START_REF] Forest | Some links between recent gradient thermo-elasto-plasticity theories and the thermomechanics of generalized continua[END_REF], [START_REF] Forest | Hypertemperature in thermoelastic solids[END_REF]).

In recent years, the motivation for deriving non-classical heat conduction models in the mechanics field stemmed from trying to understand the presence of thermal fluctuation fields in heterogeneous materials with a microstructure. Specifically, the authors in [START_REF] Forest | Hypertemperature in thermoelastic solids[END_REF] postulate the existence of a free energy function that has an added dependence on the gradient of the entropy density variable. Based on this enhanced free energy, an enhanced heat equation was derived containing a term with a characteristic length related to material parameters that can account for scale-size thermal effects in micro-heterogeneous bodies. Finally, all of the above theoretical or computational non-classical approaches seem to have found some validation in recent experimental work where evidence of size-dependent thermal effects were reported in heterogeneous materials (see, e.g., [START_REF] Fehér | Size effects and beyond-Fourier heat conduction in room-temperature experiments[END_REF], [START_REF] Fehér | On the evaluation of non-Fourier effects in heat pulse experiments[END_REF]).

A c c e p t e d V e r s i o n

In this work, we commence with a higher-gradient heat equation model. By introducing physically relevant scaling arguments related to the absolute size of the constituents, we point out a new length scale parameter that models scale-dependent thermal effects (see, e.g., [START_REF] Fehér | Size effects and beyond-Fourier heat conduction in room-temperature experiments[END_REF][START_REF] Fehér | On the evaluation of non-Fourier effects in heat pulse experiments[END_REF]). Hence, by arguing as in the work of [START_REF] Nika | Effective medium theory for second-gradient nonlinear elasticity with chirality[END_REF], the solution of the enriched microscopic problem can be seen as a vanishing-viscosity solution that coincides with the classical homogenized solution of (1.3). However, unlike in the classical case, the local solutions w k satisfy a higher-gradient local problem and, hence, possess better regularity properties than classical local solutions. What is remarkable is, that this aforementioned higher regularity of the local solution compensated by the mild assumption that the data are of the form f := div F F F , F F F ∈ L 3 (Ω, R 3 ) allows us to prove Proposition 1.1 under the minimal assumptions of a uniformly Lipschitz open set Ω and non-smooth material coefficients.

We have organized the paper as follows: In Section 2 we explain in detail the scaling argument we employ, present the enriched microscopic model, provide some motivation for its use, and prove some general qualitative results as they pertain to existence and uniqueness of solution as well as the variational nature of the problem. In Section 3 we state the main results, we discuss their consequences, and demonstrate symmetry relations for the higher-gradient effective coefficients as well as explore their variational structure. Section 4 is dedicated to proving the main results in Section 3. Finally, we reserve Section 5 for some discussion and remarks.

2 Problem set-up

Scale-dependent thermal effects

In this paragraph, we postulate the modified heat flux in (1.11) and derive an additional length scale parameter that encapsulates the size-dependent thermal effects in the context of an idealized periodic microstructure. The approach is motivated from generalized continuum mechanic theories (see, e.g., [START_REF] Forest | Milieux continus généralisés et matériaux hétérogènes[END_REF]) since it is well understood that the effective properties of heterogeneous materials can depend not only on the volume fraction of the phases or their geometrical distribution but also on the absolute size of the constituents (see, e.g., [START_REF] Forest | Asymptotic analysis of heterogeneous Cosserat media[END_REF]).

We assume that our working domain Ω has a periodic microstructure with period ℓ and overall characteristic length L. This introduces a natural periodicity of the tensors K and M. We scale all the parameters in the model, including the material parameters, the following way:

x x x * = x x x L , u * (x x x * )= u(x x x) U , (2.1)
where U := max z z z∈Ω |u(z z z)|. While for the material parameters we can define the normalized tensors,

K * =K/K, M * =M/M, (2.2)
A c c e p t e d

V e r s i o n

where

K:= max z z z∈Y ℓ |K(z z z)|, M:= max z z z∈Y ℓ |M(z z z)| with Y ℓ :=(-ℓ/2, ℓ/2]
3 the periodic cell characterizing the body Ω. We can now introduce an additional length scale relation between K and M as follows:

M = ℓ 2 TE K. (2.
3)

The characteristic intrinsic length ℓ TE accounts for scale-size effects in heat conduction problems of heterogeneous media. When it is identically zero Fourier's classical law of heat conduction is recovered.

Remark 2.1. We chose to work with the above scaling because of its physical interpretation. The scaling indicates that the size of the heterogeneities are comparable to the order of the period. Naturally, one could consider a different scaling than the one proposed above and arrive at different mascroscopic heat conduction problems. We will not address other type of scalings here. We leave the interested reader to consult the work in [START_REF] Nika | Effective medium theory for second-gradient nonlinear elasticity with chirality[END_REF] for different type of scalings and their influence on effective equations within the context of second gradient elasticity. Moreover, the works in [START_REF] Aifantis | Internal length gradient (ILG) material mechanics across scales and disciplines[END_REF], [START_REF] Forest | Some links between recent gradient thermo-elasto-plasticity theories and the thermomechanics of generalized continua[END_REF] provide a theoretical framework as to why such internal lengths are needed and how they can account for thermal scale-size effects, while the works in [START_REF] Fehér | On the evaluation of non-Fourier effects in heat pulse experiments[END_REF], [START_REF] Fehér | Size effects and beyond-Fourier heat conduction in room-temperature experiments[END_REF] provide an experimental justification.

Thus, the scaled heat flux becomes:

q i = 3 j=1   K U L K * ij ∂u * ∂x * j + 3 k,l=1 KU ℓ TE L 2 ∂ ∂x * j M * ijkl ∂ 2 u * ∂x * k ∂x * l   . (2.4)
If we use the notation, q q q * :=(KU) -1 q q q then we have a scaled form of the heat flux,

q q q * = 1 L K * ∇ * u * + ℓ TE L 2 div * (M * :∇ * ∇ * u * ) . (2.5)
To simplify our presentation, henceforth, we absorb the factor 1/L into the notation K * . Moreover, we remark that since the coefficients K and M are Y ℓ periodic the corresponding scaled coefficients K * and M * are Y * periodic where Y * := ℓ L Y with Y :=(-1/2, 1/2] 3 . Finally, hereon, if no confusion arises we will drop the * notation in order to expedite our presentation.

The microscopic problem

We consider a material with a periodic microstructure of period ε:=ℓ/L≪1 occupying a region Ω ⊂ R 3 . The region Ω that the heterogeneous material occupies is assumed to be a uniformly Lipschitz open set (see [START_REF] Demengel | Functional Spaces for the Theory of Elliptic Partial Differential Equations[END_REF]Definition 2.65]). The exterior boundary component will be denoted by Σ:=∂Ω while the vector n n n will denote the unit normal on Σ pointing in the outward direction. The ε periodic problem, generated by defining the nondimensional number ε as the ratio of ℓ/L, will permit us to obtain an effective equation when ε → 0. However, unlike in classical homogenization problems, different cases ought

A c c e p t e d V e r s i o n

to be considered depending on how the intrinsic length scale ℓ TE scales with ℓ (or L).

Here, since we are interested in recovering Fourier's classical law of heat conduction as an effective limit, we will only consider the scaling,

ℓ TE /ℓ ∼ 1.
(2.6)

The physical meaning of the above scaling, is that the intrinsic length ℓ TE is comparable with the length of the heterogeneities. Naturally, other type of scalings are possible, however, we will not address other cases here. We refer the reader to [START_REF] Nika | Effective medium theory for second-gradient nonlinear elasticity with chirality[END_REF] for different type of scalings in the context of generalized continuum mechanics.

Therefore, under the scaling in (2.6), the (generalized) heat flux becomes,

q q q ε = K( x x x ε ).∇u ε + ε 2 div M( x x x ε ):∇∇u ε . (2.7)
The microscopic problem is then characterized by the following equation and boundary conditions,

-div K( x x x ε ).∇u ε -ε 2 div M( x x x ε ):∇∇u ε = f in Ω, ε 2 M( x x x ε )∇∇u ε :n n n ⊗ n n n = 0 on Σ, u ε = 0 on Σ, (2.8)
where f is some given source that belongs in L 2 (Ω). We remark, that prescribing a homogeneous Dirichlet boundary condition, as is usually the case, is no longer sufficient. We require, additionally, to prescribe a zero heat flux for what we refer to as a normal double heat flux that is directly related to the spatial retardation coefficient M. It describes that the internal structure of the material has an effect on the heat flow. It is connected to the choice of the scaling ε 2 in front of the higher gradient term in the flux. It plays here the role of a natural boundary condition as one can see at the level of the weak formulation of the problem. It is also a convenient choice of boundary condition as it cancels any surface or line integrals. In direct analogy with second-gradient elasticity, such a term is referred to as normal double traction (see, in particular, [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, I: Théorie du second gradient[END_REF], [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF], [START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF]).

Notation and assumptions

-We employ the Einstein notation of repeated indices unless otherwise stated.

-Throughout the work we assume that the uniform strong ellipticity condition holds, i.e., there exist positive (generic) constants c 1 and c 2 such that: c 1 |w w w| 2 ≤ w w w.K(x x x).w w w ≤ c 2 |w w w| 2 , c 1 |w w w| 2 |q q q| 2 ≤ w w w ⊗ q q q:M(x x x):w w w ⊗ q q q ≤ c 2 |w w w| 2 |q q q| 2 , (2.9) for all w w w, q q q ∈ R 3 -{0 0 0}.

A c c e p t e d V e r s i o n

Auxiliary formulas

For the readers' convenience and for the expediency of our results, we introduce certain formulas that we will make use of in obtaining the variational formulation of (2.8). These formulas, among others, can also be found in [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, I: Théorie du second gradient[END_REF]Appendix].

For any smooth enough scalar function ξ defined on Σ or on a neighborhood of Σ, the tangential and normal components of ∇ξ are,

(∇ξ) τ = -n n n × (n n n × ∇ξ)=∇ξ -(∇ξ) n n n n, (∇ξ) n :=∇ξ.n n n. (2.10)
Moreover, we introduce the surface gradient of ξ using the projection operator Π:=I-n n n⊗n n n,

∇ s ξ=(I -n n n ⊗ n n n)∇ξ=Π∇ξ. (2.11)
Thus, we can write down a useful integration by parts formula on surfaces,

Σ ∇ s ξ ds = Σ ξ(div n n n)n n n ds + ∂Σ ξν ν ν dℓ, (2.12)
where

ν i = ϵ ϵ ϵ ijk t j n k (2.13)
are the components of the unit normal vector on ∂Σ and tangent to Σ, t t t is the unit tangent vector to ∂Σ, and ϵ ϵ ϵ is the Levi-Civita tensor. Lastly, we remark, the jump term on (2.12) is on a ridge, i.e. the line on Σ where the tangent plane of Σ is discontinuous.

Remark 2.2. The above formulas are used with a high degree of frequency in emulsions and capillary fluids (see e.g. [START_REF] Nika | Rate of convergence for a multiscale model of dilute emulsions with non-uniform surface tension[END_REF]). We refer the reader to the appendix of reference [START_REF] Germain | La méthode des puissances virtuelles en mécanique des milieux continus, I: Théorie du second gradient[END_REF] or [START_REF] Germain | The method of virtual power in continuum mechanics. Part 2: Microstructure[END_REF] for an excellent exposition of the above formulae and related topics.

Weak formulation

The primary setting for for the variational formulation of (2.8) is the space H 2 (Ω) ∩ H 1 0 (Ω) where the Sobolev space H 2 (Ω) is a Hilbert space with norm,

∥u∥ H 2 (Ω) = ∥u∥ 2 L 2 (Ω) + ∥∇u∥ 2 L 2 (Ω,R 3 ) + ∥∇∇u∥ 2 L 2 (Ω,R 3×3 ) 1/2 . (2.14)
Since the variational formulation of (2.8) is not a standard one, we write down the details for the readers convenience using the notation introduced in Section 2.2.2. Hence, if we multiply (2.8) by a test function v ∈ {ϕ ∈ C ∞ (Ω) | ϕ = 0 on Σ} and integrate by parts several times (including integration by parts on surfaces using formula (2.12)) we obtain,

A c c e p t e d V e r s i o n

- Σ K ij ( x x x ε )∂ x j u ε -ε 2 ∂ x k M ikpq ( x x x ε )∂ 2 xpxq u ε n i v ds + Ω K ij ( x x x ε )∂ x j u ε ∂ x j v dx x x-ε 2 Σ M ikpq ( x x x ε )∂ 2 xpxq u ε n k n m n i ∂ xm v dσ + ∂Σ M ikpq ( x x x ε )∂ 2 xpxq u ε n k Π iℓ ν ℓ v dl- Σ Π ℓm ∂ xm (M ikpq ( x x x ε )∂ 2 xpxq u ε n k Π iℓ )v dσ +ε 2 Ω M ikpq ( x x x ε )∂ 2 xpxq u ε ∂ 2 x i x k v dx x x= Ω f v dx x x.
(2.15)

Using the fact that we have imposed a homogeneous Dirichlet boundary condition and a zero normal double heat flux for the spatial retardation on Σ, we can see that the variational formulation (in vectorial form) reduces to the following: Find

u ε ∈ H 2 (Ω) ∩ H 1 0 (Ω) such that, Ω K( x x x ε )∇u ε .∇v dx x x+ε 2 Ω M( x x x ε )∇∇u ε :∇∇v dx x x= Ω f v dx x x, (2.16)
for all v ∈ H 2 (Ω) ∩ H 1 0 (Ω). Remark 2.3. The weak form in (2.15) also provides a way for recovering the strong form of problem (2.8) in the sense of distributions. For more details, the interested reader can consult [START_REF] Nika | Effective medium theory for second-gradient nonlinear elasticity with chirality[END_REF].

Remark 2.4. Existence and uniqueness of a solution in (2.16) that belongs in H 2 (Ω) ∩ H 1 0 (Ω) is a matter of applying the Lax-Milgram lemma together with Poincaré's inequality in H 1 0 (Ω) and the assumptions regarding the ellipticity of the tensors in Section 2.2.1. Hence, immediately, we can derive the following estimate from (2.16),

∥u ε ∥ 2 H 1 (Ω) + ε 2 ∥∇∇u ε ∥ 2 L 2 (Ω,R 3×3 ) 1/2 ≤ c(Ω) ∥f ∥ L 2 (Ω) .
(2.17)

Variational formulation

The weak solution to (2.16) can be classified as the unique minimum of the functional J ε (θ),

u ε = arg min θ∈H 2 (Ω)∩H 1 0 (Ω) J ε (θ), (2.18)
where

J ε (θ):= 1 2 Ω K( x x x ε )∇θ.∇θ dx x x+ 1 2 Ω ε 2 M( x x x ε )∇∇θ:∇∇θ dx x x- Ω f θ dx x x. ( 2

.19)

A standard computation of the variational derivative of J ε will recover (2.16) and the Euler-Lagrange equations in succession.

A c c e p t e d V e r s i o n

3 Main results

Homogenization via unfolding Γ-convergence

We define the following domain decompositions (see [START_REF] Ciorȃnescu | Éclatement périodique et homogénéisation[END_REF][START_REF] Damlamian | An elementary introduction to periodic unfolding[END_REF][START_REF] Ciorȃnescu | The periodic unfolding method in homogenization[END_REF][START_REF] Ciorȃnescu | The Periodic Unfolding Method. Theory and Applications to Partial Differential Problems[END_REF]): Regarding properties of the unfolding operator, the reader can consult [START_REF] Ciorȃnescu | Éclatement périodique et homogénéisation[END_REF][START_REF] Damlamian | An elementary introduction to periodic unfolding[END_REF][START_REF] Ciorȃnescu | The periodic unfolding method in homogenization[END_REF][START_REF] Ciorȃnescu | The Periodic Unfolding Method. Theory and Applications to Partial Differential Problems[END_REF].

K - ε := ℓ ∈ Z d | ε(ℓ + Y ) ⊂ Ω , Ω - ε := int ∪ ℓ∈K - ε ε(ℓ + Y ) , Λ - ε := Ω\Ω - ε . (3.1) ε ε Ω Ω - ε 0 ∈ R 3 ε x ε ∈ Z 3 εy ∈ εY x ∈ Ω
Definition 3.2. ([8, Def. 12 and Prop. 14, pg. 458]) Let F ε : L p (Ω) → R be a sequence of functionals and

F : L p (Ω × Y ) → R (p > 1)
. We say that F ε unfolding Γ-converges to F if for all u ∈ L p (Ω × Y ):

1. For every sequence u ε ∈ L p (Ω) such that T ε (u ε ) ⇀ u in L p (Ω × Y ) one has, F (u) ≤ lim inf ε→0 F ε (u ε ); (3.4)
A c c e p t e d V e r s i o n

2. there exists a sequence

u ε ∈ L p (Ω) such that T ε (u ε ) ⇀ u in L p (Ω × Y ) and
F (u) = lim ε→0 F ε (u ε ). (3.5)
Theorem 3.1. The sequence of functionals J ε unfolding Γ-converge in the weak topology of H 2 (Ω) ∩ H 1 0 (Ω) to J eff as ε → 0,

J ε Γ → J eff , (3.6)
where ) is what one would expect as a result from the homogenization of second order linear elliptic equations while the second part is new, specific to the inclusion of the second-gradient thermal effects.

J eff (θ):= 1 2 Ω K eff ∇θ.∇θ dx x x - Ω f θ dx x x, (3.7
Theorem 3.2. Let u ε and u be the unique minimizers of (2. [START_REF] Fehér | Size effects and beyond-Fourier heat conduction in room-temperature experiments[END_REF]) and (3.7), respectively. Moreover, for ũ1 ∈ L s (Ω), let f :=div

F F F , F F F ∈ L 3 (Ω, R 3 ) then, ∇u ε -∇u -∇ y û(•, • ε ) → 0 in L 2 (Ω, R 3 ), (3.10) as ε → 0 with û(x x x, y y y) = -w k (y y y) ∂u ∂x k (x x x) + ũ1 (x x x). ( 3

A c c e p t e d V e r s i o n

3.2 Remarks on the homogenized coefficients

Symmetry

Since we work within a variational framework, the homogenized coefficients inherit the symmetry that is imposed on them from the framework. If one were to obtain the effective tensor through a multiple scale expansion then the tensor would have the following, nonsymmetric, form:

K eff ik := Y K ij (y y y)(δ jk - ∂w k ∂y j
) dy y y.

(3.12)

Naturally, the two forms, (3.8) and (3.12), are equivalent (since the tensors K and M are assumed to be isotropic). The proof is virtually identical to the classical case (see [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF]), however, we feel it should be included here since it is a higher-gradient generalization of the classical case:

Multiply with a test function v ∈ H 2 per (Y ) equation (3.9) and integrate by parts to obtain,

0 = Y M mjpq ∂ 2 w k ∂y p ∂y q ∂ 2 v ∂y m ∂y j dy y y + Y K mj ∂(w k -y k ) ∂y j ∂v ∂y m dy y y. (3.13)
Selecting v = w i we have,

0 = Y M mjpq ∂ 2 w k ∂y p ∂y q ∂ 2 w i ∂y m ∂y j dy y y + Y K mj ∂(w k -y k ) ∂y j ∂w i ∂y m dy y y. (3.14)
If we add and subtract coordinate y i on the w i term of the second integral we obtain,

Y K ij (y y y)(δ jk - ∂w k ∂y j ) dy y y = Y M mjpq ∂ 2 w k ∂y p ∂y q ∂ 2 w i ∂y m ∂y j dy y y + Y K mj ∂(w k -y k ) ∂y j ∂(w i -y i )
∂y m dy y y.

(3.15)

Variational characterization of the higher-order effective tensor

From equation (3.8), we can see that K eff is symmetric as well, and is determined in its entirety from the knowledge of the quadratic form K eff ξ ξ ξ • ξ ξ ξ for any constant vector ξ ξ ξ ∈ R

-div y (K (∇ y w ξ -ξ ξ ξ) -div y (M:∇ y ∇ y w ξ )) =0 in Y, y y y → w ξ (y y y) is Y -periodic. (3.17)
Using completely standard techniques stemming from the calculus of variations we can write for our case,

K eff ξ ξ ξ • ξ ξ ξ= inf v∈H 2 per (Y ) Y K(y y y)(∇ y v-ξ ξ ξ).(∇ y v-ξ ξ ξ) dy y y+ Y M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y y , (3.18) with H 2 per (Y ):={v ∈ H 2 (Y ) | v is Y -periodic}.
It follows immediately that the effective coefficients are elliptic (they satisfy the Legendre condition), i.e., K eff ξ ξ ξ • ξ ξ ξ ≥ K|ξ ξ ξ| 2 where K= Y K(y y y) dy y y.

4 Proof of the main results

Proof of Theorem 3.1

Theorem 3.1. The sequence of functionals J ε unfolding Γ-converge in the weak topology of H 2 (Ω) ∩ H 1 0 (Ω) to J eff as ε → 0,

J ε Γ → J eff , (3.6)
where

J eff (θ):= 1 2 Ω K eff ∇θ.∇θ dx x x - Ω f θ dx x x, (3.7)
K eff ik := Y K(y y y)(∇ y w k -e e e k
).(∇ y w i -e e e i ) + M(y y y)∇ y ∇ y w k :∇ y ∇ y w i dy y y,

with

w k ∈ H 2 (Y ), k = 1, 2, 3 the Y -periodic local solution to, -div y (K(y y y)(∇ y w k -e e e k ) -div y (M(y y y)∇ y ∇ y w k )) =0 in Y. (3.9)
Proof. Define,

E ε (θ):= 1 2 a(θ, θ) if θ ∈ H 2 (Ω) ∩ H 1 0 (Ω) +∞ if θ ∈ L 2 (Ω)\H 2 (Ω) ∩ H 1 0 (Ω), (4.1)
with A c c e p t e d V e r s i o n

a(θ, θ):= Ω K( x x x ε )∇θ.∇θ dx x x+ Ω ε 2 M( x x x ε )∇∇θ:∇∇θ dx x x, (4.2)
and

E eff (θ):= 1 2 Ω K eff ∇θ.∇θ dx x x if θ ∈ H 1 0 (Ω) +∞ if θ ∈ L 2 (Ω)\H 1 0 (Ω). (4.3)
It suffices to show that E ε Γ → E eff where the energy E ε is finite.

Γ-limit inferior inequality.

Let u ε ∈ H 2 (Ω) ∩ H 1 0 (Ω) be a solution to (2.8) then based on the estimates in (2.17) and properties of the periodic unfolding operator T ε (see Definition 3.1) we have,

• T ε (u ε ) ⇀ u in L 2 (Ω, H 2 (Y )) • T ε (∇u ε ) ⇀ ∇ x u + ∇ y û in L 2 (Ω, H 1 (Y, R 3 )) • T ε (ε∇∇u ε ) ⇀ ∇ y ∇ y û in L 2 (Ω × Y, R 3×3 ). Then, for all φ φ φ ∈ L 2 (Ω, C 1 per (Y, R 3 )) we have, 0≤ 1 2 Ω K( x x x ε )(∇u ε -φ φ φ(x x x, x x x ε )).(∇u ε -φ φ φ(x x x, x x x ε )) dx x x + 1 2 Ω ε 2 M( x x x ε )(∇∇u ε -∇φ φ φ(x x x, x x x ε
)):(∇∇u ε -∇φ φ φ(x x x, x x x ε )) dx x x.

(4.4)

Opening up the parentheses in the above expression, we obtain Since the space L 2 (Ω, C 1 per (Y, R 3 )) is dense in L 2 (Ω × Y, R 3 ), the above inequality holds for a sequence of regular functions of the form φ φ φ ε (x x x, y y y) Y -periodic in y y y with the following convergence properties, A c c e p t e d V e r s i o n

E ε (u ε )≥ Ω K( x x x ε )∇u ε .φ φ φ dx x x- 1 2 Ω K( x x x ε )φ φ φ.φ φ φ dx x x + Ω ε 2 M( x x x ε )∇∇u ε :∇φ φ φ dx x x- 1 2 Ω ε 2 M( x x x ε )∇φ φ φ:∇φ φ φ dx x x = Ω×Y K(y y y)T ε (∇u ε ).T ε (φ φ φ) dy y ydx x x - 1 2 Ω×Y K(y y y)T ε (φ φ φ).T ε (φ φ φ)
• T ε (φ φ φ ε ) ⇀ ∇ x u + ∇ y û in L 2 (Ω × Y, R 3 ) • T ε (ε∇φ φ φ ε ) ⇀ ∇ y ∇ y û in L 2 (Ω × Y, R 3×3 ).
Upon extracting a (non-relabeled) sub-sequence, we obtain in the limit inferior of (4.5),

lim inf ε→0 E ε (u ε )≥ 1 2 Ω×Y K(y y y)(∇ x u+∇ y û).(∇ x u+∇ y û)
+M(y y y)∇ y ∇ y û:∇ y ∇ y û dy y ydx x x

≥ 1 2 Ω inf v∈H 2 per (Y ) Y K(y y y)(∇ x u+∇ y v).(∇ x u+∇ y v) +M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y ydx x x = 1 2 Ω K eff ∇ x u.∇ x u dx x x
Γ-limit superior inequality.

We will construct the recovery sequence for smooth functions initially and then use a diagonalization argument to complete the proof. Define the sequence,

Let u u u ∈ H 2 (Ω) ∩ H 1 0 (Ω). Relying on the density of C ∞ 0 (Ω) into H 2 (Ω) ∩ H 1 0 (Ω)
u ε (x x x) = u(x x x) + εû( x x x ε ). (4.8)
Then,

• T ε (u ε ) → u in L 2 (Ω, H 2 (Y )) • T ε (∇u ε ) → ∇ x u + ∇ y û in L 2 (Ω, H 1 (Y, R 3 )) • T ε (ε∇∇u ε ) → ∇ y ∇ y û in L 2 (Ω × Y ).

A c c e p t e d V e r s i o n

Thus,

E ε (u ε ) = 1 2 Ω K( x x x ε )∇u ε .∇u ε +ε 2 M( x x x ε
)∇∇u ε :∇∇u ε dy y ydx x x = 1 2 Ω×Y K(y y y)T (∇u ε ).T (∇u ε )+ε 2 M(y y y)T (∇∇u ε ):T (∇∇u ε ) dy y ydx x x.

(4.9) Hence, passing to the limit as ε → 0 in the expression above we obtain,

lim ε→0 E ε (u ε ) = Ω×Y (K(y y y)(∇ x u + ∇ y û).(∇ x u + ∇ y û)+M(y y y)∇ y ∇ y û):∇ y ∇ y û) dy y ydx x x = 1 2 Ω inf v∈H 2 per (Y ) Y K(y y y)(∇ x u+∇ y v).(∇ x u+∇ y v)+M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y ydx x x = 1 2 Ω K eff ∇ x u.∇ x u dx x x (4.10)
We can conclude the proof using a density and diagonalization argument. for all ϕ ∈ H 1 0 (Ω). In the sense of distributions we can recover,

-div K eff ∇u = f in D(Ω), u = 0 on ∂Ω, (4.12)
with K eff given in (3.8).

The above equation is precisely the limit problem that has been obtained countless times in the homogenization literature. The difference in this work is that the local solution used to construct the effective tensor K eff satisfies a higher-gradient problem with w k ∈ H 2 per (Y ) (and not only in H 1 (Y )). This newly found local solution regularity is combined with the theorem below to prove the classical corrector convergence result.

The next theorem is included only for completion purposes and is a consequence of G. Stampacchia's interpolation theorem [START_REF] Stampacchia | The spaces L (p,λ) , N (p,λ) and interpolation[END_REF]. Practically, it allows one to control the L p norm of the gradient of the weak solution of a second order elliptic system with constant coefficients, for p ∈ [2, ∞), by controlling the integrability of the data. We refer the A c c e p t e d V e r s i o n interested reader to [4, Thm. 3.28, pg. 76] (and references therein) for a modern treatment on the subject. Our aim here is to apply the theorem below to the second order elliptic system with constant coefficients in (4.12) to obtain the necessary integrability for the homogenized solution to compensate for the remaining regularity that is needed in the classical corrector convergence theorem.

-div(A∇u) = -div F F F , (4.13)
where the constant coefficients A ijαβ satisfy the Legendre-Hadamard condition, and

F F F ∈ L p (Ω, R 3 ) for some 2 ≤ p < ∞. Then, ∇u ∈ L p (Ω, R 3×3 ) and ∥∇u∥ L p (Ω,R 3×3 ) ≤ c ∥F F F ∥ L p (Ω,R 3 ) ( 4
(Ω), let f :=div F F F , F F F ∈ L 3 (Ω, R 3 ) then, ∇u ε -∇u -∇ y û(•, • ε ) → 0 in L 2 (Ω, R 3 ), (3.10)
as ε → 0 with û(x x x, y y y) = -w k (y y y) ∂u ∂x k (x x x) + ũ1 (x x x).

Proof. The main idea of the proof essentially amounts to a compensated regularity argument, where the burden of regularity is split (un-equally) between the local solution and the homogenized solution. We argue as follows: by the assumption on the body force f we can apply Theorem 4.1 and obtain ∇u ∈ L 3 (Ω, R 3 ). Moreover, since the local solution belongs in H 2 per (Y ), w k ∈ H 2 per (Y ), we extend it by periodicity to the entire space. By applying standard Sobolev embedding theory we can obtain ∇ y w k ∈ L 6 (R 3 , R 3 ) (see e.g. [START_REF] Demengel | Functional Spaces for the Theory of Elliptic Partial Differential Equations[END_REF]Theorem 2.31]). Therefore, the above compensated regularity argument makes the expression ∇u ε -∇u -∇ y û(•, • ε ) well defined in the L 2 norm and hence, we can prove Proposition 9.12 in [10, pg. 185] directly. In what follows, we provide the actual steps of the argument that is given in [10, Proposition 9.12, pg. 185] with the appropriate adjustments that the second-gradient thermal effects introduce.

Embarking from the ellipticity of the second order tensor K we are led to, where the last inequality is due to the ellipticity condition of the tensor M. Hence, opening up the expression above, using the structure of equation (2.8), and the symmetry of the tensors K and M we have, By unfolding, we can see that we are able to pass to the limit in each expression due to estimate (2.17 (4.17)

On the last step we made use of equation (4.10).

Discussion

By enriching a microscopic model with higher-gradient contributions, we have obtained an upscaled Fourier's law of heat conduction, which exhibits thermal effects at varying length scales through the effective conductivity tensor. The inherent higher gradient regularity of the enriched microscopic model is transferred to the local problem, enabling us to use a compensated regularity argument and establish the classical corrector result in homogenization with less stringent regularity assumptions. By and large, we expect that the results can be extended to the time-dependent case, benefiting from the natural A c c e p t e d V e r s i o n gradient flow structure of the problem (see [START_REF] Nika | A gradient system for a higher-gradient generalization of Fourier's law of heat conduction[END_REF]). Most likely the techniques of evolutionary Γ-convergence (see e.g. [START_REF] Mielke | On evolutionary Γ-convergence[END_REF]) would be applicable. Naturally, the technique here can be extended to perforated domains by adjusting the unfolding operator as in [START_REF] Ciorȃnescu | The periodic unfolding method in domains with holes[END_REF]. The challenge is to handle interface terms.

Figure 3 . 1 : 2 )

 312 Figure 3.1: Unfolding operator on a periodic grid. Let [z z z] Y = (⌊z 1 ⌋, ⌊z 2 ⌋, ⌊z 3 ⌋) denote the integer part of z z z ∈ R 3 and denote by {z z z} Y the difference z z z-[z z z] Y which belongs to Y . Regarding our multiscale problem that depends on a small length parameter ε > 0, we can decompose any x x x ∈ R 3 using the maps [•] Y : R 3 → Z 3 and {•} Y : R 3 → Y the following way (see Fig. 3.1 (right)),x x x = ε x x x ε Y + x x x ε Y . (3.2) Definition 3.1. ([13, Def. 2.1, pg. 1588]) For any Lebesgue measurable function φ on Ω we define the periodic unfolding operator by,

. 9 ) 3 . 1 .

 931 ) K eff ik := Y K(y y y)(∇ y w k -e e e k ).(∇ y w i -e e e i ) + M(y y y)∇ y ∇ y w k :∇ y ∇ y w i dy y y, (3.8)with w k ∈ H 2 (Y ), k = 1, 2, 3 the Y -periodic local solution to, -div y (K(y y y)(∇ y w k -e e e k ) -div y (M(y y y)∇ y ∇ y w k )) =0 in Y. (3Remark The first part of (3.8

3 .

 3 Using definition (3.8) one can check that, K eff ξ ξ ξ • ξ ξ ξ = Y K(y y y)(-ξ ξ ξ+∇ y w ξ ).(-ξ ξ ξ+∇ y w ξ ) dy y y + Y M(y y y)∇ y ∇ y w ξ :∇ y ∇ y w ξ dy y y, (3.16) where w ξ is the solution of the following cell problem, A c c e p t e d V e r s i o n

dy y ydx x x + Ω×Y ε 2 M

 2 (y y y)T ε (∇∇u ε ):T ε (∇φ φ φ) dy y ydx x x-1 2 Ω×Y ε 2 M(y y y)T ε (∇φ φ φ):T ε (∇φ φ φ) dy y ydx x x ε→0 → Ω×Y K(y y y)(∇ x u + ∇ y û).φ φ φ dy y ydx x x-1 2 Ω×Y K(y y y)φ φ φ.φ φ φ dy y ydx x x + Ω×Y M(y y y)∇ y ∇ y û:∇ y φ φ φ dy y ydx x x-1 2 Ω×Y M(y y y)∇ y φ φ φ:∇ y φ φ φ dy y ydx x x.

  we can suppose without loss of generality that u ∈ C ∞ 0 (Ω). Furthermore, let û ∈ C 2 per (Y ) be a minimizer of, inf v∈H 2 per (Y ) Y K(y y y)(∇ x u+∇ y v).(∇ x u+∇ y v) dy y y+M(y y y)∇ y ∇ y v:∇ y ∇ y v . (4.6) Existence of such a minimizer is shown using classical arguments of coercivity and lower semi-continuity of, v → inf v∈H 2 per (Y ) Y K(y y y)(∇ y v-ξ ξ ξ).(∇ y v-ξ ξ ξ) dy y y+ Y M(y y y)∇ y ∇ y v:∇ y ∇ y v dy y y . (4.7)

Remark 4 . 1 .

 41 The Euler-Lagrange equation of (3.7) is, Ω K eff ∇u.∇ϕ dx x x = Ω f ϕ dx x x (4.11)

Theorem 4 . 1 .

 41 Let u ∈ H 1 0 (Ω, R 3 ) be a weak solution of the Dirichlet problem

A

  c c e p t e d V e r s i o nc 1 ∇u ε (x x x)-∇u(x x x) -∇ y û(x x x, ε (x x x)-∇u(x x x)-∇ y û(x x x, x x x ε )).(∇u ε (x x x)-∇u(x x x)-∇ y û(x x x, ε (x x x) -∇u(x x x) -∇ y û(x x x, x x x ε )).(∇u ε (x x x) -∇u(x x x) -∇ y û(x x x, ε -∇ y ∇ y û(x x x, x x x ε)):(ε∇∇u ε -∇ y ∇ y û(x x x,

c 1 2 L 2 (

 122 ∇u ε (x x x)-∇u(x x x)-∇ y û(x x x, x x x ε ) y û).(∇u+∇ y û)dx x x ε :∇ y ∇ y ûdx x x+ Ω M( x x x ε )∇ y ∇ y û:∇ y ∇ y ûdx x x.

). Thus, lim ε→0 c 1

 1 ∥∇u ε -∇u-∇ y û∥ 2 L 2 (Ω,R 3 ) ≤ Ω×Y f udy y ydx x x -Ω×Y K(y y y)(∇u+∇ y û).(∇u+∇ y û)dy y ydx x x -Ω×Y M(y y y)∇ y ∇ y û:∇ y ∇ y ûdy y ydx x x=0.

  .11) Remark 3.2. We note in the above theorem, the assumptions that ∇ y w k ∈ L r and ∇u ∈ L s with 1/r + 1/s = 1/2 are no longer needed. Furthermore, if we compare (3.11) with G.

	Griso's result in (1.9), we see that (3.11) is more accessible to computations, hence, more
	practical. Note that û ∈ L 2 (Ω, H 1 per (Y )).

  .14) A proof of this theorem can be found in [26, Thm. 7.1, pg. 138] and [4, Thm. 3.29, pg. 79 and discussion on Sect. 3.5, pg. 78].

4.2 Proof of Theorem 3.2 Theorem 3.2. Let u ε and u be the unique minimizers of (2.19) and (3.7), respectively. Moreover, for ũ1 ∈ L s

Interestingly, unlike the case of Fickian diffusion, a rigorous microscopic interpretation of the classical linear Fourier law (and of its variations) done at the level of stochastically interacting particles is still lacking in spite of sustained efforts in both probability and statistical physics communities (cf. e.g.[START_REF] Bonetto | Fourier's law: a challenge for theorists[END_REF]). It would be interesting to see, for instance, whether the simplified mechanistic approach from[START_REF] Ruelle | A mechanical model for Fourier's law of heat conduction[END_REF] unveiling connections between Hamiltonian chains of coupled mechanical systems and the linear Fourier law can be eventually extended to capture at least some of the second gradient effects postulated in (1.10).

Acknowledgements

The authors gratefully acknowledge the financial support by the Knowledge Foundation (project nr. KK 2020-0152). Moreover, we would like to express our gratitude to the anonymous reviewers for their many comments, suggestions, and corrections. They have, without a doubt, improved the quality of the manuscript.


Luciano Drozda
email: drozda@cerfacs.fr

Denis Matignon
email: denis.matignon@isae-supaero.fr

Diffusive series representation for the Crandall model of acoustic impedance

Keywords: Cylindrical tubes, Acoustic impedance, Bessel functions, Diffusive representation Mathematics Subject Classification MSC 35Q35, MSC 33C10

): 2321-2330, 2016) in their study of pulse propagation within blood vessels and by Colombaro et al. (Meccanica 52(4-5):825-832, 2017) in their work on the Bessel models of linear viscoelasticity.

Introduction

Porous media are well-known for their sound absorption properties and can be found in soundproofing systems of recording studios as well as in the walls of aircraft turbojets, for example. Modeling their physical behavior has been the subject of many papers based either on Biot's theory of poroelasticity or on the Crandall model of acoustic impedance for cylindrical tubes. The latter was notably used by Maa [START_REF] Maa | Potential of microperforated panel absorber[END_REF] for the study of micro-perforated panel (MPP) absorbers, devices that became popular for their robustness in harsh environments such as aircraft engine nacelles.

The expression of the Crandall model of impedance contains a ratio of (transcendental) Bessel functions. Recently, MATLAB codes have been developed to approximate Bessel functions with complex arguments. However, the historical difficulty in accurately evaluating them near zeros [START_REF] Hart | Computer approximations[END_REF] has motivated the acoustic community to make use of polynomial approximations of the model at low and high frequencies instead. In an intermediate range of frequencies, polynomial ad hoc functions whose limiting values respect these approximations are typically used. Alternative expressions have been derived based on continued fraction approximations [START_REF] Bilbao | Passive models of viscothermal wave propagation in acoustic tubes[END_REF][START_REF] Thompson | Analog model for thermoviscous propagation in a cylindrical tube[END_REF]. Such an approach proves useful for time-domain numerical simulations but might present either slow convergence to the exact form (i.e., high approximating order is required if wideband accuracy is desired) [START_REF] Thompson | Analog model for thermoviscous propagation in a cylindrical tube[END_REF] or lack of flexibility: indeed, a reoptimization has to be performed at each change of temperature, tube radius, frequency range of interest or cost function [START_REF] Bilbao | Passive models of viscothermal wave propagation in acoustic tubes[END_REF].

In this paper, an exact representation of the impedance of an acoustic tube as a series of high-pass filters in the frequency domain is proposed. The method applied herein relies on the computation of a timedomain kernel h(t) from its Laplace transform H(s) by means of the Bromwich integral:

where c is any real abscissa such that all singularities of H lie in the left-half plane ℂ - 0 ∶= {s ∈ ℂ, ℜe(s) < 0} . It was also used by Giusti and Mainardi [START_REF] Giusti | A dynamic viscoelastic analogy for fluid-filled elastic tubes[END_REF] in their study on the effects of blood viscosity on the propagation of a pulse within a vessel. They showed that the time-domain counterpart of a so-called relaxation memory function could be represented as a Dirichlet series depending only on the zeros of the Bessel function of the first kind of order 0. This technique was also applied by Colombaro et al. [START_REF] Colombaro | A class of linear viscoelastic models based on Bessel functions[END_REF] in their study of the Bessel models of linear viscoelasticity. More mathematical discussions around it can be found in other works by these authors [START_REF] Colombaro | Bessel models of linear viscoelasticity[END_REF][START_REF] Giusti | On infinite series concerning zeros of Bessel functions of the first kind[END_REF].

In the scope of the acoustics of MPPs, such an alternative representation is of great numerical interest, since the evaluation of the ratio of Bessel functions is no longer necessary and suitable approximations of the impedance on a frequency range of interest can be easily found by series truncation.

The paper is structured as follows: Sect. 2 is dedicated to briefly recall the assumptions behind the Crandall model and the derivation of the original impedance

h(t) = 1 2j lim R→∞ ∫ c+jR c-jR
e st H(s) ds , expression in the Laplace domain; Sect. 3 shows how to achieve a series representation for it; in Sect. 4 the authors compare the original expression and its (truncated) series representation over a wide frequency range and also perform a parametric analysis on the latter. In the Appendix, a plot of the principal value of the argument of the non-dimensional model locating its zeros and poles is provided.

Physical modeling

In this section, the original model is written in the Laplace domain, along with a recall on its assumptions.

A useful result

Irving B. Crandall [START_REF] Crandall | Theory of vibrating systems and sound[END_REF] derived an expression for the transfer impedance of a cylindrical perforation in which an acoustic wave is travelling and undergoes viscous effects at the wall (see Fig. 1). It can be stated as follows, ∀s ∈ ℂ�P Cr :

where (s) ≡ √ s 0 ∕ is a complex-valued wavenum- ber related to viscous diffusion; z 0 is the air character- istic impedance, 0 the air density, its dynamic viscosity, is the distance between the two sections used for defining the transfer impedance (any effect due to the finite length of the perforation is neglected), d its diameter, and I 1 and I 0 are the modified Bessel func- tions of the first kind of order 1 and 0, respectively. In Eq. ( 1), P Cr corresponds to the set of poles of ẐCr .

(1)

ẐCr (s) = z 0 2 (s) 1 - 2 (s) d∕2 I 1 ( (s) d∕2) I 0 ( (s) d∕2) -1
, Fig. 1 A cylindrical hole with cross-section S. The unit vector points outward the fluid medium (and inside the hole)

Derivation of the result

This section is dedicated to show how Eq. ( 1) arises from the solution of the linearized momentum equation for an axisymmetric Stokes flow. In a Stokes regime, fluid flow is incompressible, such that the density 0 is assumed to be a constant and the mass continuity Navier-Stokes equation reads where U is the flow velocity field. Moreover, in a Stokes flow, the inertial forces are small when compared to the viscous forces. By eliminating the inertial terms from the momentum Navier-Stokes equation, one obtains where P is the pressure field. Under linear acoustics, one can express velocity and pressure fields respectively as U ≡ + u and P ≡ P 0 + p , where ||u|| ≪ || || and p ≪ P 0 are small perturbations to the uniform velocity and to the uniform pressure P 0 of a steady medium. This leads to the following linearized mass and momentum equations:

As stated in Sect. 2.1, the Crandall model of impedance concerns cylindrical perforations, such that the use of a cylindrical coordinate system is preferred. Besides, an axisymmetric flow parallel to the cylinder walls is considered, such that the acoustical velocity field is independent from the azymuthal direction and can be expressed as u ∶= u(r, z, t)e z , where r is the radial distance from the center of the cylindrical tube, z is the height and t is the time coordinate. Additionally, the acoustical pressure field is assumed to be uniform over the cylinder cross-section and must vary linearly along its length. Under these conditions, Eqs. ( 4)-( 5) read ( 2)

⋅ U = 0 , (3) 0 t U = -P + ∇ 2 U , (4) ⋅ u = 0 , (5) 0 t u = -p + ∇ 2 u . (6) z u = 0 ,
where (p 2p 1 ) is the (time-dependent) pressure dif- ference between both sides of the tube. It is worth noting that the viscosity is normally independent of this pressure difference. Liquids under extreme pressure can experience an increase in viscosity, but such a situation is outside the scope of the model studied here. Finally, a no-slip boundary condition at the cylinder walls is considered.

By applying the Laplace transform, the solution may be expressed as, for r ∈ [0, d∕2], z ∈ [0, ]: from which Eq. ( 1) can be deduced, since where ⟨⋅⟩ S is the cross-section surface average operator.

Indeed, since where identity ( * ) stems from relations §9.6.26 and §9.6.27 in Abramowitz and Stegun [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF], one obtains from which Eq. ( 1) proves straightforward. By making use of the recurrence relation (see §9.6.26 in Abramowitz and Stegun [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF]) [START_REF] Giusti | A dynamic viscoelastic analogy for fluid-filled elastic tubes[END_REF]

0 t u = -(p 2 -p 1 )∕ + ∇ 2 u , (8) û(r, z, s) = - 1 2 (s) 1 - I 0 ( (s) r) I 0 ( (s) d∕2) ( p2 -p1 ) , (9) ẐCr ∶= - 1 z 0 ( p2 -p1 ) ∕ ⟨û⟩ S , (10)
⟨I 0 ( (s) r)⟩ S ∶= 1 S � S I 0 ( (s) r) dΩ , = 1 d 2 ∕4 � d∕2 0 I 0 ( (s) r) 2 r dr , = 1 2 (s) d 2 ∕8 � (s) d∕2 0 s � I 0 (s � ) ds � , ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ ( * ) ≡[s � I 1 (s � )] (s) d∕2 0 (11) ⟨û⟩ S = - 1 2 (s) � 1 - 2 (s) d∕2 I 1 ( (s) d∕2) I 0 ( (s) d∕2) � , ( 12
) I n-1 (s) -I n+1 (s) = 2n s I n (s) ∀s ∈ ℂ, ∀n ≥ 1,
it is possible to rewrite Eq. ( 1) more simply as

Series representation

In this section, the derivation of an alternative expression of the Crandall impedance as a series of high-pass filters depending only on the zeros of the Bessel function J 2 is presented. In Sect. 3.1, for the sake of sim- plicity, one derives the series representation for a nondimensional formulation of the resulting expression. Finally, Sect. 3.2 presents the series representation of the Crandall impedance in both Laplace and time domains after re-injection of the physical constants into the non-dimensional formulation.

Derivation for a non-dimensional model

For the sake of simplicity, the derivation is first performed on a model which replaces all physical constants in ẐCr by unity:

The following statements hold:

Statement 1.
All the poles of Ẑnon-dim are negative reals (i.e., P non-dim ⊂ ℝ -). Let ∈ P non-dim . Rela- tion §9.6.3 by Abramowitz and Stegun [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF] ensures that 0 = I 2 ( √ ) = -J 2 (j √ ) . Let J ∶= {(j 2,k ) k∈ℕ } be the set of increasing positive zeros of J 2 (note that {(±j 2,k ) k∈ℕ } ∪ {0} is the set of all zeros of J 2 ). By choosing the principal branch of the square root function (i.e., ∀s ∈ ℂ, Arg √ s ∈ ] -∕2, ∕2] ), one concludes that Arg j √ ∈ ]0, ] and so j √ ∈ J . Therefore, Statement 2. For large modulus, the following asymptotic expansion holds:

Ẑnon-dim (s) ∼ s + 2 √ s + o(1)
, hence lim s→∞ Ẑnon-dim (s) = +∞ . From the asymp- totic expansion for large arguments of I n presented in (13)

ẐCr (s) = z 0 2 (s) I 0 ( (s) d∕2) I 2 ( (s) d∕2) ∀s ∈ ℂ�P Cr . ( 14
) Ẑnon-dim (s) ∶= s I 0 ( √ s) I 2 ( √ s) (s ∈ ℂ�P non-dim ). ( 15
) P non-dim = {(-j 2 2,k ) k∈ℕ } ⊂ ℝ -. §9.7.
1 in Abramowitz and Stegun [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF], it is possible to conclude that

In complement to Statement 1., it is worth noting that Ẑnon-dim is a single-valued function (and so it does not present any branch cut), in spite of its dependence on the square root function. This can be seen from the ascending series representations of I 0 ( √ s) and I 2 ( √ s) . Indeed, ( 16)

non-dim ∼ s 1 + 1 8 √ s + o 1 √ s 1 -15 8 √ s + o 1 √ s , ∼ s 1 + 1 8 √ s + o 1 √ s 1 + 15 8 √ s + o 1 √ s , ∼ s 1 + 2 √ s + o 1 s , ∼ s + 2 √ s + o(1) s→∞ ----→∞ .
Fig. 2 Closed integration path typically used for the Laplaceinversion of meromorphic functions with no branch cut. The integral over Γ R vanishes when R → ∞ in case Jordan's lemma applies (see §1.4 in Duffy [START_REF] Duffy | Transform methods for solving partial differential equations[END_REF]). The singularities of the function are highlighted by crosses following relation §9.6.10 in Abramowitz and Stegun [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF]. On the one hand, Statement 1. ensures that Ẑnon-dim can be Laplace-inverted. On the other hand, Statement 2. makes it impossible to apply Jordan's lemma on the Bromwich integral (see Fig. 2). One is led rather to invert Ẑmod (s) ∶= Ẑnon-dim (s)

s -1 ≡ I 0 ( √ s) I 2 ( √ s) -1 . Since Ẑmod (s) ∼ 1∕ √ s as s → ∞
, this function is eligible to the application of Jordan's lemma.

The set of poles of Ẑmod is P non-dim ∪ {0}.

For t > 0 , the computation of the Bromwich inte- gral by Cauchy's residue theorem then leads to as the associated time-domain kernel.

Alternatively, in the Laplace-domain:

for s ∈ ℂ�P non-dim .

The original expression [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] involves the square root function s ↦ √ s , which usually gives rise to a branch cut on ℝ -. But, thanks to (19), the singu- larities of Ẑnon-dim are reduced to a countable set of poles localized in ℝ -.

Back to the original model

By introducing the physical constants c 1 ∶=

𝜌 0 z 0 > 0 and c 2 ∶= √ 𝜌 0 𝜂 d 2 > 0 , one shows that for s ∈ ℂ�P Cr , with P Cr = {(-j 2 2,k ∕c 2 2 ) k∈ℕ } ⊂ ℝ -, ẐCr admits the representation (17) Ẑnon-dim (s) ∶= s I 0 ( √ s) I 2 ( √ s) = 2 ∑ ∞ k=0 (s∕4) k (k!) 2 ∑ ∞ k=0 (s∕4) k k!(k+1)! , (18)
z mod (t) = Res( Ẑmod (s) e st ;s 0 = 0)

+ ∞ ∑ k=1 Res( Ẑmod (s) e st ;s k = -j 2 2,k ) , = 8 + 4 ∞ ∑ k=1 e -j 2 2,k t , (19) Ẑnon-dim (s) = s + 8 + 4 ∞ ∑ k=1 s s + j 2 2,k ,
In the time domain, with input u, one has:

where the causal impulse response reads:

Numerical application

Series truncation for the numerical use of representation ( 20) is investigated: first, in Sect. 4.1 the effect of truncation is analyzed, and second in Sect. 4.2, a parametric study is provided when the model parameters vary in a physically meaningful range of interest.

Comparison between the series representation and its approximation

By taking a finite number of zeros of the Bessel function J 2 and with the parameter values indicated in Table 1, Bode diagrams of ẐCr and of its approxi- mated representation through series truncation are illustrated in Fig. 3a. Besides, the real and imaginary parts of ẐCr (j ) and of its approximated form are plotted in Fig. 3b. Finally, Fig. 3c shows a plot of the absorption coefficient taken as 1 -| (j )| 2 , with (j ) ∶= ẐCr (j )-1 ẐCr (j )+1

. The curves related to the exact representation ẐCr are obtained from Eq. (20) using MATLAB with its full numerical precision. On the one hand, at low frequencies ( f < 100 Hz), it is possible to see that the approximated form (20)

ẐCr (s) = c 1 c 2 2 Ẑnon-dim (c 2 2 s) = c 1 s + 8 c 1 c 2 2 + 4 c 1 c 2 2 ∞ ∑ k=1 s s + j 2 2,k ∕c 2 2 . ( 21
) (z Cr ⋆ u)(t) = c 1 u � (t) + 8 c 1 c 2 2 u(t) + 4 c 1 c 2 2 ∫ t 0 h(t -𝜏) u � (𝜏) d𝜏, (22) h(t) ∶= ∞ ∑ k=1 exp - j 2 2,k c 2 2 t , for t > 0 .
approaches the exact one even when any series term is considered, which is expected since the series is a sum of high-pass filters. On the other hand, at high frequencies ( f > 100 Hz), discrepancies are observed for a small number of filters ( k = {0, 1} ), specially in the plots of the real part (Fig. 3b) and of the absorption coefficient (Fig. 3c). It is however worth noting that only k = 10 filters already provide a satisfactory approximation, which enforces the low-numericalcost aspect of such a method. Indeed, this approach does not require any optimization procedure as compared to other works [START_REF] Bilbao | Passive models of viscothermal wave propagation in acoustic tubes[END_REF][START_REF] Hélie | Representations with poles and cuts for the time-domain simulation of fractional systems and irrational transfer functions[END_REF][START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF].

Asymptotic behaviors for the Crandall impedance can also be extracted from the representation (20): at very low frequencies ( f ≪ 10 Hz), ẐCr behaves like

s ↦ 8c 1 ∕c 2 2 , whereas s ↦ c 1 s + (2c 1 ∕c 2 ) √ s + 8c 1 ∕c 2 2
approximates it at frequencies beyond 10 4 Hz. Nev- ertheless, interest on its asymptotic behavior at very high frequencies is likely to be diminished by the fact that the incompressibility condition underlying the Crandall model is violated at such a wide range of frequencies.

The Stokes number St ∶= �c 2 √ s� = c 2 √ (a nondimensional quantity proportional to the ratio of the tube diameter to the acoustic boundary layer thickness) varies from 0.305 to 9.65 in the frequency range considered for the numerical experiments ( f ∈ [10 1 , 10 4 ] Hz -3 decades). For Stokes num- bers above value St crit = 10 , viscothermal losses could be neglected and the flow would be considered as adiabatic, a situation out of the scope of the model's assumptions.

Parametric study

Finally, a parametric analysis on the truncated series representation of ẐCr for varying values of perforation length ( ), diameter (d) and temperature has been performed. More precisely, by taking the parameter values from Table 1 as reference, Bode diagrams, real/imaginary part plots and absorption coefficient It was noticed that the model is very robust to changes in the perforation length . This is expected as the impedance ẐCr varies linearly with this param- eter. Varying the temperature in the reported range also has small impact in the plots. A more complex dependence is observed for d, which motivates Fig. 4, showing the relative errors on the Bode plots of the approximated form with respect to the exact one for different numbers of series terms ( k = {1, 5, 9}).

For a given value of k, the observed trend in the magnitude plots can be explained as follows: the Table 1 Parameter values for numerical simulations comparing expressions (1) and (20) for the Crandall acoustic impedance of a cylindrical tube filled with air; 0 , and z 0 are physical constants for air at 20 larger the diameter d, the shorter the stopband of the model's transfer function (see Fig. 5), and since the truncated series better approximates the exact form in the stopband (low frequencies) than in the passband (high frequencies), the error on magnitude must overall increase with the perforation diameter.

For the phase plots, however, the same trend is no longer observed, to which clear explanation proves difficult to find. Additionally, as k increases, the errors on magnitude and phase decrease for all values of d, as expected.

Conclusions

In this paper, a new expression for the acoustic impedance of a cylindrical tube undertaking an axisymmetric and normally incident unsteady Stokes flow is presented. The exact representation of the Crandall impedance as a series of high-pass filters depending only on the zeros of the Bessel function J 2 is of practical interest, since passage to time domain for numerical simulations is made easier than with the classical expression, and such a development can be embedded into the well-known Maa model for MPPs used in sound absorption mechanisms. Moreover, in the frequency domain, the use for numerical simulations of ad hoc connection functions linking the polynomial asymptotic behaviors of the Crandall impedance at high and low frequencies is no longer required.

Finally, the method applied herein relates to retrieving the so-called oscillatory-diffusive representations of causal kernels extensively studied in Hélie and Matignon [START_REF] Hélie | Representations with poles and cuts for the time-domain simulation of fractional systems and irrational transfer functions[END_REF], in Monteghetti et al. [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] and in Moufid et al. [START_REF] Moufid | Energy analysis and discretization of the time-domain equivalent fluid model for wave propagation in rigid porous media[END_REF] that can be used for other models of acoustics materials as well, such as the equivalent fluid ones for porous media (e.g., the Johnson-Champoux-Allard-Pride-Lafarge model [START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF]). (see §9.6.26 and §9.6.27 in Abramowitz and Stegun [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF]), it is easy to notice the following: between two given consecutive zeros of I 0 ( √ s) , a sign change of I 1 ( √ s) takes place thanks to Eq. ( 24) and to the continuity of Bessel functions; as I 2 ( √ s) and I 1 ( √ s) assume opposite signs at a root of I 0 ( √ s) because of Eq. ( 23), one concludes that I 2 ( √ s) undertakes a sign change in between the two consecutive zeros of I 0 ( √ s) and hence a root of I 2 √ s) (i.e., a pole of Ẑnon-dim ) must lie in between.

Poles-zeros interlacing seems to be a common feature for diffusive kernels with positive weights. This can be seen on the most simple rational ones: for instance, for some 𝜉 1 > 0 , 𝜉 2 > 0 , 0 ≤ ≤ 1 , the following convex combination of low-pass filters where ∶= (1 -) 1 + 2 , presents poles at s = {-1 , -2 } and a zero in between at s = -.

ĥ ∶ ℂ�{-1 , -2 } → ℂ s ↦ s + 1 + 1 - s + 2 ≡ s + (s + 1 )(s + 2 ) ,
An example of similar behavior among irrational functions is given by H (s) ∶= tanh( √ s∕ )∕ √ s for some 𝜀 > 0 , discussed in Mignot et al. [13]. Moreo- ver, this paper shows a very interesting fact: the function's intertwined poles and zeros in ℝ -become a branch cut in the limiting process leading towards 0, the limit transfer function being H 0 (s) ∶= 1 √ s .

Fig. 3

 3 Fig. 3 In (a): Bode diagrams in the frequency range [10 1 , 10 4 ] Hz (3 decades) for the exact and approximated representations of the transfer function related to the Crandall model of impedance; k = {0, 1, 10} zeros of the Bessel function J 2 have been used for the computation of the truncated series sum. In (b): Plots of the real and imaginary parts of ẐCr and its truncated series approximated representations. In (c): Absorption coefficient as a function of the frequency ▸

⋅ m - 3 1.82 ⋅ 10 - 5 1 Fig. 4

 310514 Fig.[START_REF] Colombaro | A class of linear viscoelastic models based on Bessel functions[END_REF] Plots of the relative errors on magnitude (top) and phase (bottom) of the approximated series representation with respect to the exact form for different numbers of series terms ( k = {1, 5, 9})

Fig. 5

 5 Fig. 5 Magnitude plots highlighting the evolution of the stopband of ẐCr as the diameter d increases, for k = 9 series terms

Acknowledgements

The authors would like to thank Estelle Piot from ONERA (the French Aerospace Lab) and Florian Monteghetti from ISAE-SUPAERO for the fruitful discussions.

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix: Poles and zeros in the non-dimensional model

In this Appendix, we are interested in the poles-zeros location in the complex plane of the non-dimensional model Ẑnon-dim .

Figure 6 shows a contour plot of the principal value of the argument of Ẑnon-dim (s) within a given domain of the complex plane. The set of poles is P non-dim ∶= {(-j 2 2,k ) k∈ℕ } ⊂ ℝ -while the set of zeros is {(-j 2 0,k ) k∈ℕ } ⊂ ℝ -. This plot suggests that the poles and zeros of Ẑnon-dim are intertwined. Indeed, using the following relation Publisher's Note Springer Nature neutral with regard to jurisdictional claims in published maps and institutional affiliations.