\documentclass[12pt]{article}
\usepackage[margin=2.5cm, top=2.5cm, headsep=1.5cm]{geometry}
\usepackage{multirow}%
\usepackage{amsmath,amssymb,amsfonts}%
\usepackage{mathrsfs}%
\usepackage[title]{appendix}%
\usepackage{xcolor}%
\usepackage{textcomp}%
\usepackage{manyfoot}%
\usepackage{booktabs}%
\usepackage{algorithm}%
\usepackage{algorithmicx}%
\usepackage{algpseudocode}%
\usepackage{placeins}
\usepackage{listings}%
\usepackage{placeins}
\usepackage{soul}
\usepackage{titlesec}
\usepackage{graphicx}
\usepackage{siunitx}
\usepackage{comment}
\sisetup{
  range-phrase = \ {-}\ ,
  range-units=single 
  }
\usepackage[table]{xcolor}
\usepackage{hyperref} 
\usepackage{cleveref}
\usepackage{enumitem}
\usepackage[superscript,biblabel]{cite}
\usepackage[affil-it]{authblk}
\usepackage{times} 
\usepackage{wrapfig}
\usepackage{siunitx}
\usepackage{mathrsfs}
\usepackage{hyperref} % Load the hyperref package
\hypersetup{
    colorlinks=true,
    linkcolor=blue, 
    citecolor=blue, 
    urlcolor=blue   
}
\usepackage{indentfirst}
\usepackage{colortbl}


\usepackage[superscript,biblabel]{cite}


\providecommand{\ED}[0]{$Jm/kg$}
\providecommand{\Jms}[0]{$J/m^2$}
\providecommand{\um}[0]{$\mu m$}
\providecommand{\uN}[0]{$\mu$N}
\providecommand{\uSENB}[0]{$\mu$SENB}
\providecommand{\rp}[0]{$r_p$}
\providecommand{\Rp}[0]{$R_p$}
\providecommand{\pr}[0]{$\bar\rho$}
\providecommand{\K}[1]{$K_{#1}$}
\providecommand{\s}[1]{$\sigma_{#1}$}
\providecommand{\e}[1]{$\varepsilon_{#1}$}
\providecommand{\lc}[0]{$l_{c}$}
\providecommand{\lch}[0]{$l_{ch}$}
\providecommand{\G}[1]{$G_{#1}$}
\providecommand{\W}[1]{$W_{#1}$}
\providecommand{\J}[1]{$J_{#1}$}
\providecommand{\el}[0]{$\ell$}
\providecommand{\elb}[0]{\boldmath$\ell$\unboldmath}
\providecommand{\prb}[0]{\boldmath$\bar\rho$\unboldmath}
\providecommand{\JRc}[0]{J--R curve}
\providecommand{\Sup}[0]{Supplementary Materials}
\providecommand{\RD}[0]{rhombic dodecahedron}
\providecommand{\OT}[0]{octet-truss}
\providecommand{\IP}[0]{interpenetrating}
\providecommand{\RDL}[0]{rhombic dodecahedron lattice}
\providecommand{\OTL}[0]{octet-truss lattice}
\providecommand{\IPL}[0]{interpenetrating lattice}
\providecommand{\DC}[0]{$DC$}
\providecommand{\lowDC}[0]{Low-DC}
\providecommand{\hiDC}[0]{High-DC}
\providecommand{\V}[0]{Vero}
\providecommand{\A}[0]{Agilus30}
\providecommand{\VW}[0]{Vero White}
\providecommand{\VA}[0]{Avero}
\providecommand{\cf}[0]{$c_{f}$}
\providecommand{\keywords}[1]{\textbf{\textit{Keywords:}} #1}

\renewcommand{\thepage}{S\arabic{page}} 
\renewcommand{\thesection}{S\arabic{section}}  
\renewcommand{\thetable}{S\arabic{table}}  
\renewcommand{\thefigure}{S\arabic{figure}} 
\renewcommand{\theequation}{S\arabic{equation}} 
\setcounter{figure}{0}
\setcounter{section}{0}
\setcounter{table}{0}
\setcounter{equation}{0}
\setcounter{page}{1}

\title{\textbf{Supplementary Information} \\ \vspace{2cm} Toughness by Design: Multi-Scale Interpenetrating Lattices with Size-Enhanced Fracture Resistance}

% Interpenetrating Lattices (IPL): the battle of length scale toughness enhancement
% Toughness Amplification via Preserved Size-effects in Architected Materials
% Preserving size effects via microarchitecture
% Critical Lenth Scales for Creating Tough Architected Materials

\author[a,c,\textdagger]{Abdulaziz O. Alrashed}
\author[b,\textdagger]{Zainab S. Patel}
\author[a]{Kush Dwivedi}
\author[d]{Benjamin C. White}
\author[d]{Bryan Kaehr}
\author[d]{Brad L. Boyce}
\author[a]{Lucas R. Meza\thanks{Corresponding author: lmeza@uw.edu}} 

\affil[a]{Mechanical Engineering, University of Washington, Seattle, WA, 98195, USA.}
\affil[b]{Materials Science and Engineering, University of Washington, Seattle, WA, 98195, USA.}
\affil[c]{Mechanical Engineering Department, University of Kuwait, KW, 13060, Kuwait.}
\affil[d]{Materials, Physical, and Chemical Sciences Center, Sandia National Laboratories, USA.}
\date{}

\begin{document}

    
    \maketitle
    \titlespacing*{\section}{0pt}{0.5\baselineskip}{0.5\baselineskip}    

    \keywords{Architected Toughness $|$ Nanolattices $|$ Interpenetrating Lattices $|$ Fracture Size Effects} 
    
    \footnotetext[0]{\textsuperscript{\textdagger} Co-first authors}
    
    \pagebreak
        
    \section*{Supplementary Materials}

    
    \subsection*{Nominal strength}

        The stress intensity factor (\K{I}) can be used to obtain the nominal lattice strength following the Mode I definition:
        \begin{equation}
            K_I = \sigma_N \sqrt{\pi a} ~ f\bigg(\frac{a}{W}\bigg),
            \label{eq:stress_intensity}
        \end{equation}
        \noindent
        where \s{N} is the nominal stress and $f(a/W)$ is a geometric correction factor for a particular test condition.
        For testing done in an SENB configuration \cite{srawley1976wide}, the nominal stress is
        \begin{equation}
            \sigma_N = \frac{3 S P}{2  B  W^2},
            \label{eq:sigma_SENB}
        \end{equation}
        \noindent
        where $P$ is the applied load.
        For testing done in a CT configuration \cite{srawley1976wide}, the nominal stress is
        \begin{equation}
            \sigma_N = \frac{P}{B \sqrt{W}}.
            \label{eq:sigma_CT}
        \end{equation}
        In both cases, the stress at which yielding or failure initiates (\s{NC}) can be found by taking the peak load ($P_{peak}$) in the experiment. 
        We note that this does not take into account the effect of blunt cracks \cite{qiao2019strength}, which will invariably skew the effective strength of the lattices. 
        
        \begin{table}[ht]
            \centering
            \begin{tabular}{|c||c||c|c|c|}
                \hline
                & \textbf{Solid ($\sigma_N$)}& \textbf{RD ($\sigma_N$)}& \textbf{OT  ($\sigma_N$)}& \textbf{IPL ($\sigma_N$)}\\ \hline
                \textbf{\lowDC} & 31.7& 4.8& 8.7& 4.6\\ \hline
                \textbf{\hiDC} & 71.3& 6.6& 18.7& 12.8\\ \hline
                \textbf{\VA{}} & 3.9& 0.25& 0.6& 0.8\\ \hline
                \textbf{\V{}} & 7.9& 1.1& 1.4& 0.8\\ \hline
            \end{tabular}
            \caption{Nominal strength for macro- and nanolattices. All units are MPa.}
            \label{tab:sn}
        \end{table}



    \subsection*{Enhanced SEL Strength}
        Along with a change in ductility, there is a corresponding increase in the effective strength of a strut in an SEL framework. 
        This effective lattice strut strength (\s{N-strut}) as compared to a solid material of the same size (\s{Ns}) can be quantified using the nominal strengths in Eq. 8 as

        \begin{equation}
            \sigma_{N-strut} = \sigma_{Ns} \sqrt{\frac{1+D/r_p}{1+d/r_p}}.
            \label{eq:SEL-strengthening}
        \end{equation}

        \noindent
        From this, we posit that the strength of the struts in the macroscale lattices should increase by $\sim$5.3x for \V{} and $\sim$3.7x for \VA{}, while the strength of the \hiDC{} and \lowDC{} nanolattices would only increase by $\sim$1.8x and $\sim$1.05x, respectively.
        This increase in strength can in part be used to explain the increased \W{f} of the \V{}, \VA{} and \hiDC{} lattices, which is able to counteract the 70\% reduction in energy dissipation due to being a lower density.
        Correspondingly, because the \lowDC{} nanolattices were ductile at both the structure and the strut scale, the mechanism for energy dissipation did not change and since the solid material was already at the max strength limit, the size effect from strut size reduction and the architectural effect from IPL cannot enhance it any further and therefore, the \W{f} only scaled with density. 

        \begin{table}[ht]
            \centering
            \begin{tabular}{|c||c||c|c|c|}
                \hline
                & \textbf{Solid ($\sigma_{Ns}$)}& \textbf{RD ($\sigma_{N-strut}$)}& \textbf{OT  ($\sigma_{N-strut}$)}& \textbf{IPL ($\sigma_{N-strut}$)}\\ \hline
                \textbf{\lowDC} & 31.7& 33.2& 33.2& 33.2\\ \hline
                \textbf{\hiDC} & 71.3& 126.8& 126.8& 126.8\\ \hline
                \textbf{\VA{}} & 3.9& 13.5& 13.8& 14.7\\ \hline
                \textbf{\V{}} & 7.9& 35.1& 36.7& 41.8\\ \hline
            \end{tabular}
            \caption{Nominal strength of the lattice strut for macro- and nanolattices. All units are MPa.}
            \label{tab:sn2}
        \end{table}


    \subsection*{Enhanced Strain to Failure (Size Effect)}

        \begin{figure}
            \centering
            \includegraphics[width=1\linewidth]{S1.JPG}
            \caption{Size effect enhanced strain to failure. A) Schematic of single edge notch tension (SENT) fracture specimen. B) Key stress and strain material parameters for Quasi-brittle materials. C) Key stress and strain material parameters for strain hardening (Ductile) materials. D) Definition of the Fracture energy ($G_f$) as area under the curve beyond the elastic region. E) Illustration of the size affected strain enhanced mechanisms from damage softening contributions. F) plot of the derived relation between strain enhancement and the plastic zone to size ratio as shown in Eq. \ref{eq:S_strainfy}.}
            \label{fig:S1}
        \end{figure}

        For SENT fracture test sample with the provided dimensions shown in Figure \ref{fig:S1}A, the nominal strain is defined as:
        
        \begin{equation}
            \varepsilon_{nominal} = \frac{\Delta CMOD}{L}
            \label{eq:enom}
        \end{equation}
        
        Where $L$ is the sample length and $\Delta CMOD$ is the crack mouth opening displacement defined as \cite{BazantBook2021}:  

        \begin{equation}
           \Delta CMOD = \delta_c + \Delta elastic
            \label{eq:cmod}
        \end{equation}

        Where $\delta_c$ is the cohesive zone displacement and $\Delta elastic$ is deflection in the remaining part of the beam.
        Therefore, taking Eq. \ref{eq:cmod} with Eq. \ref{eq:enom}, we obtain: 

        \begin{equation}
            \varepsilon_{nominal} = \frac{\delta_c + \Delta elastic}{L}
            \label{eq:enom2}
        \end{equation}

        The Cohesive zone displacement and elastic deflection are defined as \cite{BazantBook2021}:  

        \begin{equation}
            \delta_c  = \frac{2 \cdot G_f}{\sigma_u}
            \label{eq:delta}
        \end{equation}
        \begin{equation}
            \Delta elastic  = \frac{\sigma_f \cdot b}{E} \qquad or \qquad  \frac{P_{peak} \cdot b}{EA}
            \label{eq:delastic}
        \end{equation}
        
        Where $G_f$ is fracture energy (Figure \ref{fig:S1}D) taken from the stress-strain fracture response (Figure \ref{fig:S1}B) , $\sigma_u$  is tensile strength of material $\approx$   $\sigma_0$ maximum material fracture strength from SEL curve,  $\sigma_f$ is failure (fracture) stress,  $b$ is uncracked ligament length, $E$ is the Young's modulus, $P_{peak}$ is peak load, and $A$ is sample area. 
        Therefore, the nominal strain at fracture becomes: 

        \begin{equation}
            \varepsilon_{f } = \frac{\frac{2 \cdot G_f}{\sigma_0} + \frac{\sigma_f \cdot b}{E}}{L} = \frac{2 \cdot G_f}{\sigma_0 \cdot L} \ + \  \frac{\sigma_f \cdot b}{E \cdot L}
            \label{eq:ef}
        \end{equation}

        From the definition of the SEL curve and setting the transition size as the size of the plastic zone ($D_0 \approx r_p$): 

        \begin{equation}
            \sigma_f = \frac{\sigma_0}{\sqrt{1+D/D_0}} \qquad = \qquad  \frac{\sigma_0}{\sqrt{1+D/r_p}}
            \label{eq:SEL2}
        \end{equation}

        Then, by setting ( $b = D - a$ ) and taking Eq. \ref{eq:SEL2} into Eq. \ref{eq:ef}, the strain to fracture becomes: 


        \begin{equation}
            \varepsilon_{f} =  \frac{2 \cdot G_f}{\sigma_0 \cdot L} \ + \  \frac{\sigma_f }{E } \cdot \frac{D-a}{L} \ = \ \frac{2 \cdot G_f}{\sigma_0 \cdot L} \ + \  \frac{1}{E } \cdot \frac{\sigma_0}{\sqrt{1+D/r_p}} \cdot \frac{D-a}{L}
            \label{eq:ef2}
        \end{equation}

            For a Standard SENT specimen \cite{british2014method}, the ratio of sample dimensions can be set to ( $L = 10D$ and $a = D/2$ ):

        \begin{equation}
            \varepsilon_{f} =  \frac{2 \cdot G_f}{10 \cdot \sigma_0 \cdot D} \ + \  \frac{1}{E } \cdot \frac{\sigma_0}{\sqrt{1+D/r_p}} \cdot \frac{1}{20}
            \label{eq:ef3}
        \end{equation}

         Assuming a quasibrittle material ( $l_c = r_p$ ), the fracture energy can be written as:

        
        \begin{equation}
            G_f  = \frac{l_c \cdot \sigma_0^2}{E} \ = \ \frac{r_p \cdot \sigma_0^2}{E}
            \label{eq:Gf}
        \end{equation}

        Therefore, 
        
        \begin{equation}
            \varepsilon_{f} =  \frac{ r_p \cdot \sigma_0^2}{5 \cdot \sigma_0 \cdot D \cdot E} \ + \  \frac{1}{E } \cdot \frac{\sigma_0}{\sqrt{1+D/r_p}} \cdot \frac{1}{20}
            \label{eq:ef4}
        \end{equation}

        \begin{equation}
            \varepsilon_{f} =  \frac{1}{5} \cdot \frac{r_p}{D} \cdot \frac{\sigma_0}{E} \ + \  \frac{\sigma_0}{E } \cdot \frac{1}{\sqrt{1+D/r_p}} \cdot \frac{1}{20}
            \label{eq:ef5}
        \end{equation}
        
        Since Strain at the initiation of failure is ( $\varepsilon_{0} = \frac{\sigma_0}{E}$ ), then the ratio between the two strains becomes:

        \begin{equation}
           \frac{ \varepsilon_{f}}{\varepsilon_0} =  \frac{1}{5} \cdot \frac{r_p}{D}  \ + \  \frac{1}{20}  \cdot \frac{1}{\sqrt{1+D/r_p}} \cdot 
            \label{eq:efe0}
        \end{equation}

        Therefore, the size effect enhancement in strain to failure follows the relationship shown in Figure \ref{fig:S1}D and  \ref{fig:S1}E as: 

        \begin{equation}
            \frac{\varepsilon_{f}}{\varepsilon_0} \propto \bigg(\frac{r_p}{D}\bigg).
            \label{eq:S_strainfy}
        \end{equation}


        \begin{figure}
            \centering
            \includegraphics[width=1\linewidth]{S2.JPG}
            \caption{Size affected plastic zone expansion. A) Architectural plastic zone (\Rp{}) normalized by constituent material plastic zone size (\rp{}). B) Lattice unit cell size (\el{}) influence on the architectural plastic zone size driven by the analytical solution shown in equation (\ref{eq:S_Rp1}) for different hardening materials.}
            \label{fig:S2}
        \end{figure}        


    \subsection*{Size-affected Lattice Plastic Zone Expansion}

        The plastic zone size of a material is predicted to scale as

        \begin{equation}
            R_p \propto \bigg(\frac{K_{Ic}}{\sigma_y}\bigg)^2.
            \label{eq:S_Rp}
        \end{equation}

        \noindent
        Taking the relationships for lattice fracture toughness and strength scaling in Eqs. 3 and 2, respectively, we obtain 

        \begin{equation}
            R_p \propto \bar\rho^{~2(m-q)}\bigg(\frac{\sigma_{fs}}{\sigma_{ys}}\bigg)^2\ell.
            \label{eq:S_Rp_sigma}
        \end{equation}

        \noindent

        The yielding behavior of plastic hardening materials can generally be approximated by a Ramberg-Osgood model as

        \begin{equation}
            \frac{\varepsilon}{\varepsilon_y} = \frac{\sigma}{\sigma_y} + \bigg(\frac{\sigma}{\sigma_y}\bigg)^n,
            \label{eq:S_RamOs}
        \end{equation}

        \noindent
        where $n \geq 1$ is a hardening coefficient for the material.
        Taking this to the failure limit of a material \cite{tankasala2020crack}, we can rearrange and approximate the strain as

        \begin{equation}
            \frac{\varepsilon_f}{\varepsilon_y} \approx \bigg(\frac{\sigma_f}{\sigma_y}\bigg)^{\frac{2n}{n+1}}
            \label{eq:S_RamOsFail}
        \end{equation}

        \noindent
        Taking Eq. \ref{eq:S_RamOsFail} in combination with Eqs. \ref{eq:S_Rp_sigma} and 12, we obtain

        \begin{equation}
            R_p \propto \bar\rho^{~2(m-q)}\bigg(\frac{r_p}{d}\bigg)^{\frac{n+1}{n}}\ell.
            \label{eq:S_Rp1}
        \end{equation}

        Note here that $d$ and \el{} are related to each other based on the relative density. 
        This relation is plotted in Figure \ref{fig:S1}B for different hardening coefficient materials. 



    \subsection*{Lattice case:  $R_p \approx D$  (\VA{} RD)}

        The \VA{} RD is the only sample that lies on the borderline between stable and unstable transition out of all the tested samples, as it had an architectural plastic zone size almost equal to sample size ($R_p \approx D$). 
       It is a case worth highlighting which shows the effective transition from unstable to stable failure in the load-displacement curve driven by the ratio between plastic zone size to sample size. There are two sets of data that were gathered for this case with the same RD lattice design, strut and unit cell sizes, and base material constituent. However, due to the variation in sample preparation, including the nature of fabricating the mixture of \VA{} material consisting of brittle Vero white and ductile Angilus30 blend, there can be a slight difference in sample properties with the slight variation in the mixture ratio.
 
        The two tested samples had an ($R_p /D \approx1$) for set 1 and ($R_p /D \approx2$) for set 2 which reflect two distinct failure behaviors on the load-displacement curve as seen in Figure \ref{fig:S3}B. Set 1 shows small sudden drops in the load compared to a clearly more stable drop for set 2 despite both samples reaching full failure at a reasonably similar displacement value.
        A slight variation in sample preparation led to a larger plastic zone for set 2 that pushed it deeper into the stable region based on the SEL framework, and set 1 demonstrated a mixed behavior of stable and unstable failure as it falls at the borderline between the two regions. This case shows the strength of an SEL framework and the importance of quantifying the plastic zone size to design architectures that would fall into the stable domain and demonstrate enhanced stable fracture. 


        \begin{figure}
            \centering
            \includegraphics[width=1\linewidth]{S3.JPG}
            \caption{Characteristic failure behavior of the \VA{} RD specimens. The measured \Rp{} is reported alongside the data.}
            \label{fig:S3}
        \end{figure}



    \subsection*{Influence of Model Variables on the Work of Fracture}

        In the proposed model for lattice work of fracture found in Eq. 13 and shown in Figure  \ref{fig:S4}A, there are three key variables that dictate the contribution from lattice damage softening mechanisms and plastic hardening mechanisms toward the optimal work of fracture. The damage terms ($g_0 $ \& $\gamma_0$) account for the effective width of the fracture process zone, while the plastic term ($\alpha$) accounts for the hardening contribution from the plastic zone width ($H$). 
        The rest of the variables are set by material and size choices, which include elastic modulus ($E$), hardening coefficient ($n$), plastic zone (\rp{}) of the material, strut diameter ($d$), unit cell size ($l$), relative density ($\bar\rho$) and nominal stress ($\sigma_N$) at the given size.  
        The curve peaks around $d \approx r_p$, but the value and the behavior prior and after the peak is dictated by the three variables: $g_0$ , $\alpha$ and $\gamma_0$. 
        
        We first show the pure influence of the solid material damage term ($g_0$) which is the dimensionless energy release rate for a given material geometry shape and crack to width ratio (a/W) that can be calculated following the procedure described in \cite{MEFFORD2017961, QIAO2019127, CUSATIS20092163, 10.1115/1.2806812, SALVIATO201667} . Different work of fracture curves are plotted for different values of ($g_0$) as shown in Figure \ref{fig:S4}B, assuming no change in damage mechanism ($\gamma_0$ = $g_0$) or plastic mechanisms ($\alpha = 1$) from lattice architectures. 
        The term only changes the peak value of the curve around ($d \approx r_p$) while both the large size plateau ($d >> r_p$) and small size  plateau ($d << r_p$) converge toward the relative density line (25\%). 
        In Figure \ref{fig:S4}C, all the damage terms are set equal to one ($\gamma_0=$ $g_0=1$) to observe the influence of the plastic term ($\alpha$). It changes both the peak value and the small size plateau while the large size plateau of all plotted curves converges toward the density line. In Figure \ref{fig:S4}D, the lattice damage term ($\gamma_0$) is varied while the other two are set as constants equal to one ($\alpha = 1$ and $g_0=1$). The term only changes the large size plateau while the peak and small size plateau remain the same. 
        In essence, the peak value is controlled by $g_0$ and $\alpha$, large size plateau is controlled by $\gamma_0$ and small size plateau is controlled by $\alpha$.

        This provides insight into the unique behavior of the Vero IPL which was the only sample that showed damage-induced enhancement not seen in the IPL of the other studied materials. Architectural damage mechanisms will only reflect on the work of fracture when the strut size is larger than the material plastic zone ($d>r_p$) operating in the unstable brittle region. Once the strut size becomes equal to or less than the plastic zone size, the behavior transitions into the stable region. That is where plasticity mechanisms take over, rendering the damage mechanism from IPL architecture ineffective as seen in the other materials where both IPL and OT lattice samples demonstrate almost identical energy dissipation behavior further validating that size-induced plasticity is the driving factor with minimum influence from the architecture. Therefore, based on the strut size choice, the work of fracture can be enhanced using complex architectures such as IPL which spread the damage zone for large strut sizes, or using simple architectures such as OT that undergo size brittle-to-ductile transition (BDT) and expanded plasticity for small strut sizes . 



        \begin{figure}
            \centering
            \includegraphics[width=1\linewidth]{S4.JPG}
            \caption{Lattice work of fracture model. A) Analytical model derived to predict the size and lattice architecture's influence on the work of fracture from both damage softening and plastic hardening contributions. B)  Influence of the material damage term ($g_0$) on lattice work of fracture normalized over solid material's work of fracture, assuming no change in damage mechanism ($\gamma_0$ = $g_0$) or plastic mechanisms ($\alpha = 1$) from lattice architecture. C) Influence of the lattice plastic term ($\alpha$) that dictates the plastic hardening contribution from the plastic zone width ($H$) on lattice work of fracture normalized over solid material's work of fracture, assuming no change in damage mechanism. D) Influence of the lattice damage term ($\gamma_0$) on lattice work of fracture normalized over solid material's work of fracture, assuming no change in plastic mechanism.
            \label{fig:S4}}
        \end{figure}        

        
    \subsection*{Supplementary Videos }

        Supplementary Video 1: In-situ three-point bending fracture test on nano-scale RD, OT, and IPL lattices captured using a scanning electron microscope (SEM).
        
        Supplementary Video 2: Compact tension fracture test on macro-scale RD, OT, and IPL lattices captured using a digital camera.  










    
    \pagebreak


    \bibliographystyle{unsrt}
    \bibliography{ReferencesSupp}


\end{document}
    
    