Direct ink writing (DIW) is a three-dimensional (3D) printing technique exploited by researchers working in fields from scaffolds for energy applications to bioprinting. DIW's main strength is that it enables shaping advanced materials, if these materials can be formulated into complex fluids that meet the demands of the printing process. They must be extremely shear thinning soft solids, able to flow through narrow nozzles, recovering their structure upon deposition and retaining the predesigned 3D shape. Formulation design and rheology are critical, but these aspects can be overlooked due to the high specialization required. This work provides insight on the rheology and printability of complex yield-stress fluids through the study of linear and nonlinear behaviors using large-amplitude oscillatory shear rheology. We refine previous protocols and develop tools to understand the behaviors of formulations for DIW. We apply an existing mathematical framework to a library of carbon-based formulations for energy applications. Fourier transform analysis enables quantifying the onset and rising of higher harmonic contributions. Quantitative comparisons between different formulations are established using 3D harmonics maps, stress–strain plots, and material measures of nonlinearities [Fourier and Chebyshev coefficients, elastic moduli $(GM\u2032,\u2009GL\u2032$), and dimensionless index of nonlinearity ($S$)]. 3D Lissajous plots provide a qualitative alternative to interpretate the yielding transition. We create Ashby-type printability maps to guide formulation design and elucidate that non-printable formulations show distinctive features. This insight on yield-stress fluids for DIW is relevant to other applications and technologies: drilling fluids, gels, colloids, and foods.

## I. INTRODUCTION

Direct ink writing (DIW), a three-dimensional (3D) printing technique that relies on the formulation and design of complex (yield-stress) fluids, is an expanding multi-disciplinary research field with a growing number of applications, from two-dimensional (2D) materials such as graphene-based materials and transition metal carbides and nitrides (MXenes) for energy devices,^{1–4} and composites^{5,6} to self-healing polymers and gels for tissue engineering.^{7–10} The potential of these applications is huge, and this is evidenced by the large multi-disciplinary community embracing this technology. One of DIW's main pitfalls is the limited understanding of the complex fluids that must be carefully designed for successful printing. To advance the field and fully exploit these applications, it is crucial to achieve fundamental understanding of the complex fluids that underpin these techniques. DIW requires soft solids, which yield-stress fluids that can flow on demand through narrow nozzles (extremely shear thinning) and that must recover their original structure and retain the designed shape upon deposition in very short timescales. Despite the progress made in the field and the increasing number of researchers working in this area, there are still very few published studies that report an adequate and comprehensive rheological characterization of the printable formulations. Often only a brief discussion of the formulation's “viscosity” or a flow curve is provided, while some studies do not include any characterization at all. Those studies in which the rheology is reported quite often lack a clear and relevant analysis, and the quality of the data reported could be questioned. The specialization required for the formulation and characterization of printable materials limits the roll-out of this technology, compared to readily available off-the-shell additive manufacturing techniques such as inject printing, fused deposition modeling, or stereolithography. There is a need to improve the quality of the data^{11} reported and their analysis to deepen our fundamental understanding of such complex fluids. The quest on understanding yield-stress fluids is common to many other applied research areas, for example, drilling fluids,^{12} foods,^{13} biological fluids,^{14,15} and soft matter^{16,17} in general, which supports and strengthens the motivation and purpose of this work.

Complex fluids for DIW are formulated through many different approaches including Poloxamer^{®} gels (Pluronic^{®} F127 is widely used),^{10,18} colloidal suspensions,^{19} liquid crystals, 2D materials,^{1,20} and pH-responsive branched copolymer surfactants^{3} to name a few. Among these, Poloxamers^{®} F127 (referred to as F127 in this paper) is the additive most used because it is readily available, relatively inexpensive, and easy to formulate. It is also frequently used as a drug delivery system and in bioprinting.^{7–10,21} Graphene oxide (GO) suspensions^{22–26} are also used in DIW, but not as widely as F127 yet. GO and other 2D materials in suspension exhibit a unique and fascinating rheology, and they display a wide range of behaviors from Newtonian to printable in a relatively small concentration range, between 0.01 and 3 vol. %,^{27} with the added benefits that they act as multi-functional additive to aid the printing of other materials,^{22} and can also play a role as functional material upon thermal reduction to reduced GO (rGO).^{4,20,28} The number of researchers using these materials (F127 and GO) as formulation base in DIW is rapidly growing, and hence, it is necessary to address current gaps in the field, deepening our fundamental understanding of these soft solids.

Previous studies^{22,27,29,30} and reviews^{31} on this topic generally agree that complex fluids for DIW are yield-stress, shear-thinning fluids. The concept of “printability” and the design of yield-stress fluids are highly dependent on the context of each specific application. In bioprinting, for example, low stiffness [storage modulus, $G\u2032(\gamma ,\omega )$] and yield-stress ($\sigma y$) values are required to ensure cell viability,^{9} while other applications will require a more accurate control of morphological features that can only be achieved through stiff formulations. However, the need for careful design and better understanding of yield-stress fluids is common to all the applications that can be realized through DIW. Rheological studies have been done in continuous (flow ramps) and oscillatory shear, and the most reported rheological parameters are the storage modulus, $G\u2032(\gamma ,\omega )$ and yield stress, $\sigma y$ (note that sometimes it is referred to as the flow stress, $\sigma f$). For these yield-stress fluids, the $\sigma y$ values typically range between two orders of magnitude (∼100–1000 Pa), and the $G\u2032(\gamma ,\omega )$ from ∼1000 Pa up to 1 GPa. We found that data collected in continuous shear for “printable” formulations (or any other complex fluid with high yield-stress values) are accompanied by large uncertainties^{27} and other issues (such as shear fracture, shear banding, and other transient effects) that make the data unsuitable for a reliable and quantitative analysis. We found that as we approach the printability window, the flow behavior index (*n,* dimensionless) values approach 0, and uncertainty increases dramatically.^{27} An extended Cox–Merz rule can be used to correlate the frequency sweeps in oscillatory shear with continuous shear,^{27} and however, the uncertainty of this rule also increases within the printability window. We also found that oscillatory tests on printable formulations are reliable and reproducible (with very small uncertainties),^{22,27} and they are also extremely versatile to probe our materials reducing the issues we find in continuous shear. Large-amplitude oscillatory shear (LAOS) has been used within the DIW community; however, the way it is reported is often incomplete; for example, results might be provided without a description of how the test was done or a discussion of results. This makes very difficult to compare and rely on published data.

In previous work, we have performed structure evolution tests (combining time, frequency, and amplitude sweeps) and study their breakdown through large-amplitude oscillatory shear (LAOS) tests and their recovery through LAOS to small-amplitude oscillatory shear (SAOS) transitions.^{22,27} We concluded that we can describe and predict the printability of a formulation based on the following parameters obtained from LAOS tests: the storage modulus within the linear viscosity region (LVR), $G\u2032@LVR(\omega )$; the $\sigma f$ value at the solid-to-liquid transition [moduli crossover, $G\u2032\gamma ,\omega =G\u2033\gamma ,\omega $]; and what we defined as the flow transition index (FTI, dimensionless) that quantifies the yielding region.^{27} These parameters can be used in Ashby-type charts to create printability maps. However, this simple approach might lead to flawed conclusions if the LAOS data are not carefully analyzed. Although we believe that these “maps” can and will be a useful tool (and indeed, we include examples in Sec. III of this work), LAOS has much more to offer in DIW.^{32,33}

Here, we study the rheology of a library of carbon-based formulations using Pluronic gels F127 (Refs. 10 and 18) and GO suspensions^{22,27} as formulation base, which are then mixed with different active materials, including graphite (Gr, Table I) and multi-walled carbon nanotubes (MWCNTs, referred to as CNTs for simplicity, Table I). GO plays a dual role as a formulation base and active material. The formulations are designed with very specific concentrations (Table II) that have been chosen from a previous systematic study.^{34} The selection criteria were based on their final functional performance from electrical conductivity measurements, either as printed or after thermal reduction. However, their functional performance is not the focus of this work, and its discussion will not be included. Here, we focus on rheological studies to compare the fingerprint of the formulation base (F127 or GO) and how the addition of different active materials (Gr and CNTs) changes these fingerprints. A detailed characterization and discussion of the materials used are provided in Sec. II.

Raw material . | Particle size . | BET (m^{2}/g)
. | WCA (deg) . |
---|---|---|---|

Gr | ^{*}D_{50} = 18.3 ± 0.2 μm | 4.3 ± 0.1 | 40 ± 5 |

CNTs | D ∼ 105 ± 22 nm | 57.7 ± 0.5 | 95 ± 4 |

GO | Lateral flake size ∼ 4.7 ± 1.9 μm | ⋯ | 52 ± 1 |

Raw material . | Particle size . | BET (m^{2}/g)
. | WCA (deg) . |
---|---|---|---|

Gr | ^{*}D_{50} = 18.3 ± 0.2 μm | 4.3 ± 0.1 | 40 ± 5 |

CNTs | D ∼ 105 ± 22 nm | 57.7 ± 0.5 | 95 ± 4 |

GO | Lateral flake size ∼ 4.7 ± 1.9 μm | ⋯ | 52 ± 1 |

. | . | Active materials . | ||
---|---|---|---|---|

Formulations . | Formulation base (in distilled water) . | Gr . | CNTs . | GO . |

F127 | F127 (25 wt. %) | ⋯ | ⋯ | ⋯ |

GrF127 | F127 (25 wt. %) | 38 wt. % | ⋯ | ⋯ |

GrCNTsF127 | F127 (25 wt. %) | 38 wt. % | 0.5 wt. % | ⋯ |

GO | GO (5.5 wt. %) | ⋯ | ⋯ | ⋯ |

(1.5 vol. %) | ||||

GrGO | GO (1.8 wt. %) | 31.2 wt. % | ⋯ | ⋯ |

(1 vol. %) | ||||

CNTsGO | GO (3.2 wt. %) | ⋯ | 0.5 wt. % | ⋯ |

(1.5 vol. %) | ||||

Hybrid | F127 (25 wt. %) | 38 wt. % | 0.5 wt. % | 1 wt. % |

. | . | Active materials . | ||
---|---|---|---|---|

Formulations . | Formulation base (in distilled water) . | Gr . | CNTs . | GO . |

F127 | F127 (25 wt. %) | ⋯ | ⋯ | ⋯ |

GrF127 | F127 (25 wt. %) | 38 wt. % | ⋯ | ⋯ |

GrCNTsF127 | F127 (25 wt. %) | 38 wt. % | 0.5 wt. % | ⋯ |

GO | GO (5.5 wt. %) | ⋯ | ⋯ | ⋯ |

(1.5 vol. %) | ||||

GrGO | GO (1.8 wt. %) | 31.2 wt. % | ⋯ | ⋯ |

(1 vol. %) | ||||

CNTsGO | GO (3.2 wt. %) | ⋯ | 0.5 wt. % | ⋯ |

(1.5 vol. %) | ||||

Hybrid | F127 (25 wt. %) | 38 wt. % | 0.5 wt. % | 1 wt. % |

In this work, we use LAOS^{32,35} tests to characterize the yielding transition and nonlinear behaviors of the selected formulations designed for 3D printing of conductive architectures.^{34} We use Fourier transform (FT) analysis, and a mathematical framework to determine quantitative material measures of nonlinearities proposed by Ewoldt *et al.*^{36} This approach enables us to quantify the onset of nonlinear behavior (from the FT spectrum, in Sec. III B) and material measures of nonlinearities (Sec. III D), Fourier coefficients ($Gn\u2032,\u2009Gn\u2033$), Chebyshev coefficients ($en,\u2009vn$), elastic moduli ($GM\u2032,\u2009GL\u2032)$, and dimensionless index of nonlinearity ($S$), using higher FT-harmonics information. We compare our transient data analysis, with the correlation results provided by TRIOS, the software of the TA rheometer (Sec. III C). LAOS results are also used to analyze the yielding of printable and non-printable soft solids using a different perspective^{37} that correlates LAOS and steady shear of soft solids at the onset of the fluid rheological behavior. We complement these results with 3D (see the supplementary material, *S*ec. III E) and 2D *Lissajous*–*Bowditch* (L–B) plots^{38} (Sec. III F) that provide a visual approach to qualitatively interpretate quantitative material measures such as Fourier and Chebyshev coefficients. Section III F focuses on printability maps, or Ashby^{39} type charts based on data that can be obtained from a carefully analyzed LAOS test, that can be used as a tool to guide formulation design. We aim to put our results in a wide context, but to the best of our knowledge, there are no other published studies focused on the formulations we study here; or any other quantitative LAOS FT-rheology studies on 3D printable formulations. We have found some examples of studies on high-yield-stress fluids of the additives^{40} and actives materials used in this work;^{41} however, they are not formulated at comparable concentrations. Due to the relatively small number of formulations analyzed, establishing a link between nonlinear parameters and printability is beyond the breath of this study. The results presented here will demonstrate that LAOS and FT-rheology are powerful tools: (1) to understand the nonlinearity of complex fluids designed for DIW; (2) to quantify the boundaries of linear and nonlinear behaviors; and (3) to assess how the linear and nonlinear behavior change when using different ingredients in the formulation. These analyses are a promising path to study the behavior and performance of bespoke additives and to advance formulation development in the future.

## II. MATERIALS AND METHODS

### A. Materials

*Graphite powder (Gr).* It is purchased from Sigma-Aldrich 99% carbon basis with a particle diameter of <45 *μ*m as stated on material data sheet (CAS no. 7782-42-5). The graphite raw powders have been sieved through 100 *μ*m to aid the breakdown of aggregates and discard large particles before preparing the formulations. *Multi-walled carbon nanotubes powder (CNTs).* It is purchased from Sigma-Aldrich (multi-walled >90% carbon basis) with a diameter of 110–170 nm, length of 5–9 *μ*m, and density of 1.7 g/mL as stated on the material data sheet (CAS no. 308068-56-6). *Graphene oxide (GO)*. Graphene oxide solution (GO-4-2500) was purchased from *Graphenea* (Spain) and fully characterized. Some of the results are included in this manuscript; additional information can be found in previous work.^{34} The pH of GO suspensions was determined using a S210 pH meter (InLab Expert Pro-ISM), with values of ∼2.2 for the GO stock solution. *Pluronic F127*. Sigma-Aldrich, UK. CAS no. 9003-11-6.

### B. Powder characterization techniques

Dynamic light scattering (DLS) measurements were taken to measure the size distribution of the Gr powder (dispersed in distilled water) using a Malvern Mastersizer 3000. For GO, a small amount of sample was dispersed in distilled water and deposited between two glass slides for observation in an Olympus BX53 optical microscope. Lateral flake size was determined from individual (>200) flakes on several microscopy images using ImageJ, image analysis software. For CNTs, scanning electron microscopy (SEM) images were taken in a Hitachi S4800 SEM and used to determine their diameter with ImageJ. Gr and CNT powders were uniaxially pressed using an Specac Uniaxial press with a max load of 2 ton, to make pellets with an average diameter of 8 ± 1 mm. A concentrated GO suspension was deposited on a glass slide, evenly spread, and left to dry ensuring full coverage of the glass surface. Water contact angle (WCA) measurements were carried out using the sessile drop method in a Kruss DSA100 system. Gr and CNT powders were degassed for 6 h under vacuum prior nitrogen physisorption and Brunauer–Emmett–Teller (BET) surface area measurements in a Micromeritics 3 Flex. Elemental chemical analysis (Thermo Scientific Flash 2000 configured for %CHNS) on the GO commercial suspension provided the following composition: 45.9 wt. % C, 1.9 wt. % H, and 0 wt. % N with the rest being oxygen. This gives a C/O content of 0.88, which is higher than that of the GO used in previous works.^{22,27}

### C. Materials characterization results

Although the active materials (Gr, CNT) have similar densities (∼1.7, 2.2 mg/ml) and composition (all carbon based), their shapes, sizes, wettability, and specific surface area are diverse playing a crucial role in the interparticle, and particle/solvent interactions and therefore in the formulation process, 3D printing behavior, and functional performance. Each of the active materials used in this work has different intrinsic properties (such as particle size distribution, shape, specific surface, and wettability) that determine their behavior and printability when mixed with our two formulation base systems: F127 and GO. GO has a dual role, as the additive to aid the printing of active materials, and as an active material itself after thermal reduction during post-processing. Here, we include a summary (Table I, Fig. 1) with the results of their characterization using the techniques detailed in Sec. II B.

In this work, we use a commercial GO source unlike our previous published work.^{22,27} The commercial GO flakes have a considerable smaller lateral flake size (with lateral sizes ranging between ∼1 and 12 *μ*m for commercial GO, Fig. 1) compared to our previous work, in which we synthesized our own GO flakes (with lateral size ∼64 ± 40 *μ*m obtained from image analysis).^{34} Lateral flake size has an impact in formulation performance and printability, larger flakes facilitate network formation, but this can be compromised with smaller flakes.^{42} Gr particles have an even wider particle size distribution (Fig. 1, Table I), with particles ranging from ∼1 to 100 *μ*m (Fig. 1). A small number of larger particles and possibly aggregates are also present in the DLS analysis (Fig. 1). These larger particles are discarded in the sieving process.

### D. Formulation of carbon materials for DIW

#### 1. Pluronic^{®} F127 as a formulation base

Pluronic^{®} F127 is a triblock copolymer of poly-propylene oxide (PPO) flanked by two poly-ethylene oxide blocks (PEO). Pluronic^{®} solutions of block copolymers self-assemble to form micelles when reaching a critical micellar temperature (CMT). At sufficient concentrations, the triblock micelles associate above a critical gelation temperature (CGT) and form a lyotropic liquid crystalline (LLC) phase.^{21} The PPO block dehydrates with temperature,^{21} being hydrophobic at $T\xb0C>TCGT\xb0C$ and hydrophilic at $T\xb0C<TCGT\xb0C$, which is responsible for the thermo-reversible gelling behavior. The gelation temperature depends on the type of Pluronic^{®} and concentration.^{21} F127 is a nearly Newtonian liquid at low temperature (below CMT and CGT), which makes it very easy to add and mix other components in the formulations. F127 becomes a printable hydrogel $T\xb0C>TCGT\xb0C$ ∼ 18 $\xb0C$ at a concentration of 25 wt. % in water.^{4,18,43} Recent studies on F127 show that aqueous solutions above the CGT can be modeled using the Herschel–Buckley law,^{44} and that high temperatures well above the CGT result in Newtonian, or viscoelastic fluids depending on the F127 concentration.^{45} It has been established that at 25 wt. % and above CGT, F127 in water forms a normal micellar cubic phase (noted as $I1)$ in which the micelles crystallize in a cubic lattice.^{46} This cubic phase exhibit elastic properties^{47} providing a good base for printable formulations of other materials. A key benefit of Pluronic^{®} hydrogels is that they can be used as a “carrier” of powders with a wide range of properties, for example, ceramics (oxides and non-oxides),^{18,31,43} carbon-based materials,^{4} composites,^{5} and bio-inks^{7,9} among many others.

#### 2. F127 stock solution

A 25 wt. % Pluronic F127 stock solution was prepared by adding Pluronic F127 powder into a polytetrafluoroethylene (PTFE) jar of distilled water at low temperature (∼4 °C). Small amounts of powder were added gradually followed by a 2-min mixing cycle at 2000 rpm in a THIARE250 planetary mixer. The mixture was cooled down in between cycles to prevent heating during mixing and the subsequent water evaporation, in either a refrigerator or ice bath. The process was repeated until a smooth and transparent gel at room temperature was achieved. The stock solution was stored in the refrigerator at ∼4 °C.

#### 3. F127-active material formulations

Graphite (Gr) powder was first sieved below 102 *μ*m to remove agglomerates and discard large particles that could clog the nozzle or interfere in the LAOS experiments. The Gr powders (concentrations given in Table II) were added in several steps to the Pluronic F127 stock solution in a PTFE container that was kept cooled in a refrigerator at ∼4 °C. After a small amount of powder was added to the Pluronic stock solution, the PTFE container was sealed and mixed for 2 min in the planetary mixer at 2000 rpm. After each mixing cycle, the mixture was then cooled in a refrigerator or an ice bath. The formulations were subjected to at least two de-foaming cycles in the mixer for 2 min at 2000 rpm to remove air bubbles. The formulation, GrCNTsF127, that contains both Gr and CNTs, was prepared first by adding carbon nanotubes into Pluronic F127 in the cooled liquid state in a PTFE jar, following the same process as for GrF127 for several cycles until the mixture was homogeneous. Gr powders were then added following the same method as above.

#### 4. Suspensions of 2D colloids as formulation base

2D material suspensions are a new class of soft matter in between macromolecules and colloids; this is due to their large aspect ratio and their surface chemistry (e.g., distribution of different functional groups on the basal plane and edges of the flakes). 2D GO flakes are a few atoms thick but can have a lateral flake size up to hundreds of micrometers, while the combination of un-oxidized hydrophobic islands and oxygen functional groups (carboxy*l-COOH*, hydroxy*l-OH*) on the basal plane and edges results in an amphiphilic behavior, enabling the establishment of a combination of inter-particle interactions (π–π stacking) and hydrogen bonds with other materials.^{48} GO flakes suspensions in water form networks and gels that exhibit good printing behavior at relatively low concentrations,^{27} and can also be used as the only additive in printable formulations of oxide and non-oxide ceramics, polymers, and even steel.^{22} This is due to the unique combination of chemistry and morphology that GO has an amphiphilic (surfactant) 2D colloid that can establish networks through *π*–*π* stacking and hydrogen bonding,^{48} thus able to establish interactions with a wide range of particles with different properties. Here, we use GO as the only additive to prepare formulations for DIW (Table II) containing graphite particles (GrGO), multiwalled carbon nanotubes (CNTsGO), and a hybrid formulation prepared with a combination of active materials and additives (hybrid).

** GO stock solution (1.5 vol. %)** was prepared by gradually adding freeze-dried GO powder into the commercial GO suspension, followed by mixing in a THIARE250 planetary mixer for 2 min at 2000 rpm.

#### 5. GO-active material formulations

The graphite/graphene oxide and CNTs/graphene oxide (GrGO and CNTsGO, Table II) inks were prepared using GO suspensions with concentrations of 1 and 1.5 vol. %, respectively. When using the stock GO solution (1.5 vol. %) as a formulation base for Gr powders, the texture became lumpy and heterogeneous. The best formulation was achieved using only 1 vol. % GO (Table II). Sieved graphite powder was then added gradually in a stepwise manner followed by mixing for 1 min in the mixer and left to cooldown between mixing cycles. The formulations were also subjected to at least two de-foaming (which is a specific setting of the planetary mixer) cycles for 2 min using the THIARE250 mixer.

#### 6. Large-amplitude oscillatory shear (LAOS) measurements

Rheology data were collected on an ARES G2 strain-controlled rheometer (TA Instruments) using a 1 mm gap, a 40 mm stainless steel serrated parallel plate and a solvent trap. Amplitude sweeps were performed on the samples at strain values that ranged between 0.01% and 250% at a fixed frequency of 0.5 Hz. This range was adjusted for each specific formulation (*Materials and Methodologies* section in the supplementary material). These settings were chosen to probe the structure of the samples from small- to large-amplitude oscillatory shear (SAOS to LAOS). LAOS tests were performed in correlation (logarithmic sweep, with 20 points per decade) and transient mode (1024 points per cycle, 3 s delay time, and six half cycles). The correlation tests were repeated three times for all the samples in this work. Average values and standard deviation were calculated to validate our analysis from transient data. The transient data collection can be easily set up in the Trios, TA's software, and provides the raw data [strain, $\gamma t$ (input) and stress, $\sigma t$].

#### 7. Data signal processing and analysis

Fourier transform (FT) rheology enables quantifying nonlinearities in LAOS tests. For a sinusoidal stress input, $\gamma t=\gamma 0\u2009sin(\omega t)$, the total stress response (“whole wave approach”) can be represented by a Fourier series [Eq. (1)] including only odd harmonics due to the assumption of symmetry with respect to directionality of shear strain or shear rate.^{35,38} Within the linear viscoelastic regime, the stress response will only include the first harmonic (*n* = 1); $G1\u2032,\u2009G1\u2033$ become the real $(G\u2032)$ and imaginary part $G\u2033$ of the complex modulus $G*=G\u2032+iG\u2033$. As nonlinear material responses arise at larger strain, higher harmonic contributions appear and grow in the material's stress response. These nonlinearities can be captured using the higher harmonics coefficients $Gn\u2032$ and $Gn\u2033$, instead of just using the first harmonics information ($G\u2032$ and $G\u2033$) provided by the software in the correlation mode

The raw data signals (stress and strain vs time) were prepared before FT processing in MATLAB by trimming from the first maximum of strain to guarantee that all the waves have the same length and include two full cycles. The trimmed data signals were analyzed with the FFT function. The relative intensity for each harmonic $In/I1$ is calculated from the FT spectrum. Here, we build 3D maps showing the relative intensity of higher harmonics on the frequency space at increasing strain inputs. The projection of this 3D map on a strain-frequency plot enables to accurately determine the onset of nonlinearities.

#### 8. Chebyshev coefficients

There are alternative approaches to the “whole-wave” Fourier coefficient approach to analyze LAOS tests, such as stress decomposition,^{49} strain decomposition,^{13,33} and the sequence of physical processes (SPP) developed by Rogers *et al.*^{50–52} These approaches enable studying the yielding of soft materials and understanding different types of behaviors, such as the $G\u2033$ overshoot.^{53} The SPP approach enables to discern elastic and viscous processes in the waveforms.^{50–52} These other methods are valuable and useful; in particular, the SPP method could elucidate the differences between printable and non-printable formulations and we will use them in the future to better understand the physical processes during the yielding transition. As a first step in our LAOS analysis, we use the mathematical framework proposed by Ewoldt *et al.*^{36} based on Chebyshev polynomials of the first kind that are orthogonal and exhibit symmetry about $x=0$ (where $x=\gamma \gamma 0$).^{36} This approach enables us to carry out a physical interpretation of higher harmonics and determine quantitative measures of nonlinearities. The elastic ($en$) and viscous ($vn$) Chebyshev coefficients were calculated from the Fourier coefficients obtained in our frequency analysis of the stress data signals, shown in Eqs. (2) and (3), respectively

These coefficients, $en$ and $vn$, can be positive and negative, which has different physical implications. The third harmonic is generally used and interpreted as a measure of the material's nonlinear response.^{35,36,54} In the linear regime, $e3e1\u226a1$ and $v3v1\u226a1$, so $e1\u2192G1\u2032$ and $v1\u2192G1\u2033\omega $. A positive value for $e3$ represents intracycle strain stiffening of the elastic stress, while $e3<0$ indicates strain softening. A positive value for $v3$ corresponds to intracycle shear thickening of the viscous stress, and a negative value ($v3<0$) describes shear thinning.

Elastic moduli in the nonlinear regime, $GL\u2032$ and $GM\u2032$, are also used as a material measure of nonlinearities and can be determined from Fourier and Chebyshev coefficients, and in principle, based on their definition from the Lissajous–Bowditch (L–B) curves (see Sec. II D 9)

When $GL\u2032>GM\u2032$, the material is strain stiffening within a given cycle, and $GL\u2032<GM\u2032$ corresponds to strain softening. Continuing with the same framework by Ewoldt *et al.*,^{36} these elastic nonlinear moduli can be used to calculate a strain stiffening ratio ($S)$, dimensionless index of nonlinearity.^{36} This strain stiffening ratio becomes 0 for a linear elastic response, $S>0$, and indicates strain stiffening and $S<0$ strain softening

#### 9. Lissajous analysis

Material stress responses to LAOS can be visually interpreted using *Lissajous*–*Bowditch* (L–B) plots, in which the full-stress response is plotted parametrically against the strain (elastic representation) or shear rate (viscous representation). LAOS responses can be visualized in a 3D space with the stress [$\sigma (t)$], strain [$\gamma (t)$], and strain rate [$\gamma \u0307(t)$] as the coordinate axes (Figs. SO and S6, supplementary material).^{35,38,54} From these 3D spaces, we can plot different projections, such as stress vs strain [$\sigma t\u2009vs\u2009\gamma (t)$] denoted as the *elastic L*–*B curves*, and stress vs strain rate [$\sigma t\u2009vs\u2009\gamma \u0307(t)$] denoted as the *viscous L*–*B curves*.^{38} A material with a purely elastic response forms a straight line in the elastic L–B plot, while a purely viscous response will form a circle (in a suitably scaled plot), and a viscoelastic response an ellipse in the LVR (Fig. 2).^{38} Moderate and large viscoelastic nonlinearities can be easily identified by non-elliptical distortions in the L–B curves (Fig. 2).^{38} These curves provide a valuable visual approach to qualitatively compare different formulations in the absence of an FT analysis (Sec. III E). The distortions of L–B curves can be analyzed and quantified in different ways, such us determining the minimum ($GM\u2032$) and large ($GL\u2032$) higher harmonic moduli defined in Eqs. (4) and (5). $GM\u2032$ is the minimum-strain modulus, or tangent modulus at $\gamma =0$, and $GL\u2032$ is the large-strain modulus, or secant modulus at the maximum imposed strain $\gamma =\gamma 0$ (Ref. 36) (Fig. 2).

The raw data signals were smoothed (using the moving median function with a window of 20), to create Lissajous 3D spaces [$\sigma t,\u2009\gamma t,\u2009\gamma \u0307t$,^{35,38,54} included in the supplementary material (Figs. S5 and S6) for all the formulations in this work]. Smoothing the raw data avoids “kinks” in the curves plotted in the 3D space and L–B projections that sometimes appear due to measurement noise, which becomes more noticeable when calculating the strain rate [$\gamma \u0307t$] values from the raw data signals. 3D surfaces can also be analyzed using the SPP method that utilizes the Frenet–Serret theorem, which defines the position of a moving particle along a curvature in a coordinate system with a set of vectors.^{52,55} The projection of binormal vectors is evaluated to define transient elastic ($Gt\u2032$) and viscous moduli ($Gt\u2033$).^{17,55,56} In this work, we implement Fourier transform coupled with Chebyshev polynomials,^{36} using the material measures defined in Sec. II D 8, $GL\u2032,\u2009GM\u2032$ (Fig. 2). We also determine $GL\u2032,\u2009GM\u2032$ directly from the L–B curves for the F127 stock solution as an example (Sec. III D).

Yield-stress fluids such as drilling fluids exhibit similar *n* (flow index) values to DIW formulations (*n* = 0.099 for an invert-emulsion drilling fluids;^{38} n → 0 for printable GO suspensions^{27}) albeit the latter have larger storage modulus in the LVR ($G1\u2032$ in the range of kPa to GPa). Some of the material measures that can be determined from the LAOS and FT analysis for this type of complex fluids are sensitive to the data acquisition rate, resulting in larger slopes at discontinuities.^{38} The fact that the rate of deformations is not constant within a given cycle^{57} can lead to a *“strain softening/strain hardening paradox”*^{57} and discrepancies between different methodologies (Sec. III D).

A material measure that has almost no sensitivity to the data acquisition rate is the perfect dissipation ratio, *ϕ*, which is a metric to quantify how closely a particular LAOS response is to an idealized perfect plastic yield-stress behavior.^{38} This type of ideal behavior corresponds to microstructures like the DIW formulations in this work that are disrupted by certain yield stress (unique for each formulation), but that after disruption easily flow. *ϕ* compares the energy dissipated in a single LAOS cycle [area enclosed in a Lissajous curve, in gray in Fig. 2(e), $Ed=\u222e\sigma d\gamma $] to the energy that would be dissipated in a rigid, perfect plastic response with equivalent strain amplitude $\gamma 0$, and maximum stress $\sigma max$ [area delimited by the red square in Fig. 2(e), $Edpp=2\gamma 02\sigma max$].^{38} The energy dissipated in a LAOS cycle can also be calculated using the first-order viscous Fourier coefficient; thus, *ϕ* can be calculated using the following equation:^{38}

This measure is described as “well-behaved” for any arbitrary LAOS response because the strain amplitude and maximum stress are always well-defined and easily determined from the data. A perfect plastic behavior corresponds to $\varphi \u21921$, a purely elastic response to $\varphi \u21920$, and $\varphi \u2192\pi 4\u22480.785$ corresponds to a Newtonian behavior.^{38} This measure is used in this work to compare the behavior of different DIW formulations in Sec. III E.

## III. RESULTS AND DISCUSSION

### A. Printability of carbon-based formulations for DIW

We have previously reported a protocol to determine the printability of soft materials^{27} using GO suspensions with different concentrations that exhibit a wide range of behaviors from nearly Newtonian to printable. What we mean—and it is generally understood within the DIW community—when we label a formulation as printable is that we can reproduce the 3D shape that we originally designed based on a visual judgment of print quality. We also mean that during the printing of the part: (1) the formulation flows steadily without showing dye swell effects, clogging the nozzle, showing signs of phase separation, filter press effects, or segregation; (2) that the filament retains its shape; and (3) that the final part can support its own weight without deformation, slumping, or collapsing. It is possible to assess the printability of a formulation by empirically observing its performance during printing; however, we do not recommend this unscientific approach to characterize and classify complex fluids for DIW.

Using this preliminary and empirical assessment (summarized in Table III), we found that the GrGO formulation in this work is not printable (Table III). When using the syringe pumps, this formulation does not steadily flow through the tip of the nozzle. Instead, its microstructure is disrupted under the plunger displacement, resulting in filter-press effects and water segregation. We must admit that we did not expect this result, as we had successfully formulated different materials using GO as the only additive.^{22} We hypothesize the factors that are causing this behavior are: (1) using GO from a commercial source with relatively small lateral flake size (see Sec. II); (2) combined with graphite particles that are relatively large compared to the GO flake lateral size; and possibly, but not as relevant based on WCA measurements; and (3) that the GO elemental analysis shows a lower C/O ratio. It has been reported that the larger the lateral flake size, the lower concentration is needed to form a printable gel.^{42} Factors 1 and 2 are highly likely limiting the ability of (small) GO flakes to form a network with (large) Gr particles that can break down and rebuild with ease. However, this unexpected result offers an excellent opportunity to compare this formulation with printable ones in this study, and to explain its poor performance from a rheology perspective (in Secs. III B–III E). It has not been possible to formulate a printable formulation with both active materials (Gr and CNTs); however, we have been able to do so combining CNTs and GO in the absence of Gr particles. Active materials are added in different amounts with the purpose of increasing functional performance, in this case, electrical conductivity. Gr particles, with large average size and low specific surface area (Table I), are added in large amounts to formulate conductive pastes using F127 and GO as additives (Table II). On the other hand, CNTs, with dimensions in the nm scale, large surface area and large aspect ratio (Table I) are added in very small amounts (Table II) to improve functional performance.^{34} CNTs added in such small amounts can easily create conductive pathways and improve electrical conductivity.^{34,58,59} The intrinsic properties of the active materials and additives play a critical role in the microstructure, network formation, and rebuilt, and therefore on their printability. The failure of the GrGO formulation could be explained by a lack of network formation between GO and Gr particles, due to weak, or lack of, inter-particle interactions. In short, the small GO flakes from a commercial source are not a suitable approach to prepare printable GrGO formulations, but they can be successfully used to aid the printing of CNTs, and hybrid formulations (Table III). The hybrid performs well during printing, which suggests that the combination of F127 and GO does enable the formation of a network between the Gr particles and CNTs that can be broken and rebuilt to meet DIW demands.

### B. Fourier transform (FT) rheology analysis

The sequence of strain input waves through the amplitude sweep (as $\gamma 0$ increases) and the sequence of the material's stress responses are shown in Fig. 3. The highlighted region shows a full cycle in the amplitude sweep (from $\gamma =0$ all the way to $\gamma 0$, back to 0, and then backward to $\u2212\gamma 0$ and back to $\gamma =0$). These strain input waves are perfectly sinusoidal and do not show any distortion. We have also validated the quality of the input strain waves by checking their spectra, where we observed only one dominant frequency (corresponding to the oscillation frequency, $\omega 0$) for all the strain amplitudes. Figure 3(b) shows the material's response, a sequence of stress data signals. As the strain increases in the sweep so does the stress, however as $\gamma 0$ increases over certain values, the stress waves start showing distortions [highlighted by the arrows in Fig. 3(b)]. These distortions can be represented in the form of higher order ($n\u22653$) harmonics.^{36} Figure 3(c) shows the FT results for F127 stock solution as a representative example. This 3D map shows the intensities of higher harmonics with respect to the intensity of the first harmonic ($In/I1$, z axis) vs both the frequency (with ticks and labels for the odd harmonics in the x axis, odd multiples of $\omega 0=0.5\u2009Hz$) and $\gamma 0$ (strain amplitude) values on the y axis. The 3D maps compile all the FT spectra and provide a visual “snapshot” of the onset and growth of higher harmonics as $\gamma 0$ increases. The stress data signals and FT 3D maps for all the formulations in this work are included in the supplementary material (Figs. S1–S4). Distortions on the stress waveforms become evident for all the formulations as strain increases due to the contribution of higher harmonics. The best way to determine the precise onset of non-linearities for each higher harmonic is using the *x*–*y* projection FT fingerprints (Fig. 4), plotting the strain $\gamma 0$(%) on the *y axis*, vs the FT frequency ($n\xb7\omega 0$) on the *x axis* with labels for the $n$ harmonics. FT fingerprints for all the samples in this work are compared in Fig. 4, and the main results from this analysis are compiled in Table IV.

Formulations . | FT analysis . | LAOS results . | |||||
---|---|---|---|---|---|---|---|

Onset of NL, $\gamma NL$ (%) . | Higher harmonics, n, at $\gamma 0\u2009max$
. | $\sigma y$(Pa) (onset of NL)
. | $\gamma c$ (%) . | $\sigma f$(Pa) (bulk flow)
. | $G1\u2032$(Pa) (linear regime)
. | FTI $\sigma f/\sigma y$
. | |

F127 | 2.8 | 11th ($\gamma 0\u2009max$ = 250%) | 257 | 11.2 | 450 | 11 000 | ∼2 |

GO | 0.9 | Ninth ($\gamma 0\u2009max$ = 250%) | 52 | 54 | 270 | 9000 | ∼5 |

GrF127 | 0.04 | Ninth ($\gamma 0\u2009max$ = 50%) | 58 | 7.4 | 890 | 140 000 | ∼15 |

GrCNTsF127 | 0.045 | Ninth ($\gamma 0\u2009max\u2009$= 50%) | 55 | 7 | 840 | 135 000 | ∼15 |

GrGO^{a} | 0.11 | 11th ($\gamma 0\u2009max$ = 100%) | 463 | 80 | ⋯ | 1 250 000 | ⋯ |

CNTsGO | 0.8 | Seventh ($\gamma 0\u2009max\u2009$= 100%) | 80 | 56 | 410 | 14 000 | ∼5 |

Hybrid | 0.045 | 11th ($\gamma 0\u2009max$ = 50%) | 46 | 7 | 840 | 90 000 | ∼18 |

Formulations . | FT analysis . | LAOS results . | |||||
---|---|---|---|---|---|---|---|

Onset of NL, $\gamma NL$ (%) . | Higher harmonics, n, at $\gamma 0\u2009max$
. | $\sigma y$(Pa) (onset of NL)
. | $\gamma c$ (%) . | $\sigma f$(Pa) (bulk flow)
. | $G1\u2032$(Pa) (linear regime)
. | FTI $\sigma f/\sigma y$
. | |

F127 | 2.8 | 11th ($\gamma 0\u2009max$ = 250%) | 257 | 11.2 | 450 | 11 000 | ∼2 |

GO | 0.9 | Ninth ($\gamma 0\u2009max$ = 250%) | 52 | 54 | 270 | 9000 | ∼5 |

GrF127 | 0.04 | Ninth ($\gamma 0\u2009max$ = 50%) | 58 | 7.4 | 890 | 140 000 | ∼15 |

GrCNTsF127 | 0.045 | Ninth ($\gamma 0\u2009max\u2009$= 50%) | 55 | 7 | 840 | 135 000 | ∼15 |

GrGO^{a} | 0.11 | 11th ($\gamma 0\u2009max$ = 100%) | 463 | 80 | ⋯ | 1 250 000 | ⋯ |

CNTsGO | 0.8 | Seventh ($\gamma 0\u2009max\u2009$= 100%) | 80 | 56 | 410 | 14 000 | ∼5 |

Hybrid | 0.045 | 11th ($\gamma 0\u2009max$ = 50%) | 46 | 7 | 840 | 90 000 | ∼18 |

^{a}

Not printable. This formulation exhibits a max stress ($\sigma max\u223c2500\u2009Pa$) in the MAOS region before the solid-to-liquid transition, and therefore, a bulk flow plateau cannot be determined.

The FT analysis shows quantifiable differences between the two formulation additives used in this work, F127 and GO. The FT analysis also enables us to quantify the impact of the addition of active materials, Gr and CNTs. LAOS tests of DIW formulations often show a narrow LVR region, which makes its limits extremely subjective. The FT analysis provides a quantitative determination of the LVR limits and, therefore, of the yielding transition. In this study, only the intensities of harmonics, which are more than $1%$ of the intensity of the first harmonic ($I1$), is considered for the FT analysis, and anything below $0.01I1$ is considered as noise. The onset of nonlinearities (NL) is determined when $I3I1>Inoise\u2009\u223c0.01$, and takes place at $\gamma 0\u22482.8%$ for F127 and $\gamma 0\u22480.9%$ for GO (Table IV, Fig. 4). The addition of Gr particles to the F127 hydrogel matrix triggers an early onset of the third harmonic at $\gamma 0\u2009\u22480.04%$ (Fig. 4). However, CNTs added in small amounts to increase the electrical conductivity^{34} to this GrF127 formulation do not seem to change the FT spectra, or the onset of nonlinearities (Fig. 4 and Table IV). The impact of the CNTs addition is negligible compared with the intense effect that highly concentrated Gr particles in the F127 matrix have. When using GO as formulation additive, the FT map for GrGO (not printable, Table III) shows the onset of the third harmonic at small strain values, $\gamma 0\u22480.11%$ [Fig. 4(b), Table IV]; however, this result does not flag any distinctive feature that could be associated with printability or lack of, from these FT fingerprints. The FT maps for CNTsGO and GO are very similar [both in Fig. 4(b)]. Although the onset of the third harmonic takes place at slightly smaller $\gamma 0$ for CNTsGO (Table IV), evidencing that the addition of a small amount of CNTs to GO does not have a big impact on the GO network microstructure. The hybrid formulation containing all additives and active materials (Gr, CNTs, GO, and F127) has an excellent printing behavior (Table III), and very interesting functionality when subjected at different post-processing conditions.^{34} The FT harmonics fingerprint for this hybrid (Fig. 4) is similar to those for GrF127 and GrCNTsF127. The onset of nonlinearities takes place at $\gamma 0\u22480.045$% (Fig. 4, Table IV).

3D FT harmonics maps provide a useful visual interpretation of the onset and rising of nonlinearities, and the relative intensity of the different harmonics. The 2D FT fingerprints enable us to accurately determine the limits of the LVR, and the start of the yielding transition, which is a very important measure in DIW. The FT analysis results demonstrate that LAOS and FT rheology are a powerful tool to quantify the contribution of each active material and additive on nonlinear behaviors. Even if a direct link between such nonlinear behaviors and printability does not emerge yet from our results, we can already assert that the FT analysis will help us optimize the synthesis of bespoke additives (such as pH-responsive branched copolymer surfactants)^{3} and formulation design.

### C. LAOS analysis: First harmonic moduli and stress–strain curves from dynamic strain sweeps

The first harmonic moduli, $G1\u2032$ and $G1\u2033$, calculated from the Fourier coefficients from transient data, or those given by the rheometer's software (Fig. 5), enable us to classify the behavior of the formulations. For the F127 and GO formulations, we include a comparison between the first harmonic moduli values obtained from our transient data signals FT analyses (Fig. 5), with the average and standard deviation from three amplitude sweeps performed in correlation (calculated by *TRIOS* software). The results show an excellent agreement (Fig. 5), and for simplicity, we omit this comparison for the other formulations.

Plotting the recast of first harmonic moduli [Fig. 5(a)], and the stress vs strain data on a separate plot [Fig. 5(b)], help us to compare the yielding transition of our formulations. These stress–strain graphs are inspired by a previous study that correlates LAOS and steady shear of soft solids.^{37} Note that we do not include any steady shear experiments here because the results for our printable formulations are considerably uncertain.^{27} In the stress–strain plots, we include the following rheological parameters for the analysis: $\gamma C(%)$ is the critical strain (determined at the crossover point $G1\u2032=G1\u2033$ [Figs. 5(a) and 5(c)]; and $\gamma NL(%)$ is the strain at the onset of nonlinearities determined from the $\gamma 0(%)$ vs ($n\xb7\omega 0$) plots, *x*–*y* projections of the 3D FT harmonics maps (Fig. 4, Table IV). We show that the start of the yielding transition takes place at the onset of nonlinearities [$\gamma NL%,\u2009\sigma y(Pa)$, Table IV]; from this point, the slope of the stress vs strain plot starts to change, and plateaus around the critical strain [Figs. 5(b) and 5(d)]. We consider this “plateau” on the stress values as the “bulk flow” stress, noted as $\sigma f(Pa)$. The extent of the yielding transition is delimited between $\gamma NL%$ and $\gamma c(%)$. Here, we also provide a comparison of the $\gamma NL%$ and $\sigma y(Pa)$ values obtained through the quantitative FT approach (Table V) with the values that would be obtained in its absence by using the following approaches: (1) considering the end of the LVR when $G1\u2032<0.90GLVR\u2032$, or when $G1\u2032<0.95GLVR\u2032$ (this is the most commonly used method in DIW); and (2) when a change of slope is observed in the stress amplitude vs stress amplitude log –log plot [Figs. 5(b) and 5(d)].

Formulations . | $\gamma NL%$ . | |||
---|---|---|---|---|

$G1\u2032<0.95GLVR\u2032$ . | $G1\u2032<0.90GLVR\u2032$ . | Slope $\sigma \u2009Pa\u2009vs\u2009\gamma 0%$ plot . | $I3/I1>Inoise\u223c0.01$ . | |

F127 | 1.41 | 2.02 | 2.5 | 2.8 |

GO | 0.13 | 0.19 | 0.6 | 0.9 |

$\sigma yPa$ | ||||

$G1\u2032<0.95GLVR\u2032$ | $G1\u2032<0.90GLVR\u2032$ | Slope in $\sigma \u2009Pa\u2009vs\u2009\gamma 0%$ plot | $I3/I1>Inoise\u223c0.01$ | |

F127 | 148 | 199 | 237 | 257 |

GO | 11 | 17 | 43 | 52 |

Formulations . | $\gamma NL%$ . | |||
---|---|---|---|---|

$G1\u2032<0.95GLVR\u2032$ . | $G1\u2032<0.90GLVR\u2032$ . | Slope $\sigma \u2009Pa\u2009vs\u2009\gamma 0%$ plot . | $I3/I1>Inoise\u223c0.01$ . | |

F127 | 1.41 | 2.02 | 2.5 | 2.8 |

GO | 0.13 | 0.19 | 0.6 | 0.9 |

$\sigma yPa$ | ||||

$G1\u2032<0.95GLVR\u2032$ | $G1\u2032<0.90GLVR\u2032$ | Slope in $\sigma \u2009Pa\u2009vs\u2009\gamma 0%$ plot | $I3/I1>Inoise\u223c0.01$ | |

F127 | 148 | 199 | 237 | 257 |

GO | 11 | 17 | 43 | 52 |

F127 stock solution shows a type III behavior [weak strain overshoot, Fig. 5(a)], while the GO formulation shows a type I behavior [strain thinning, Fig. 5(c)].^{61} Type III behavior has been reported for other plastic soft materials, Carbopol, Xanthan Gum, and concentrated Ludox also display this behavior.^{53} This very recent work based on strain decomposition investigates plasticity and the yielding transition to differentiate between two modes of dissipation (solid-like and liquid-like),^{53} and the overshoot in the loss modulus is related to plasticity and a gradual yielding transition. This yielding transition is key in DIW, and it is directly related to printability. Our research and experience provide evidence demonstrating that printable formulations are complex fluids type I or type III;^{61} however, we do not yet see a direct correlation between type III and printability. Both formulation base systems, F127 and GO, are printable and exhibit similar values of $G1\u2032$ and $G1\u2033$ in the linear regime (∼10 000 Pa), and both show a strain thinning behavior in the nonlinear regime. However, there are also clear differences between them: the linear regime for GO is narrower; F127 shows a loss modulus overshoot; and they exhibit a clearly different critical strain $\gamma c(%)$ [Figs. 5(a) and 5(c)]. The stress–strain plots [Figs. 5(b) and 5(d)] show that both, F127 and GO, follow a similar trend with slightly different stress values (Table IV). Both show a clear stress plateau [Figs. 5(b) and 5(d)] with differences in their yielding region, which is very narrow for F127.

The values for $\gamma NL%$ and $\sigma yPa$ obtained using different methods (Table V) demonstrate that the most accurate and objective approach to determine the onset of nonlinear behavior (and therefore the limit of the LVR) is using the FT analysis. In the absence of an FT analysis, the change in the slope of the stress vs strain log–log plot [Figs. 5(b)], and Fig. 5(d) is the approach that leads to less errors in the determining $\gamma NL%$ and $\sigma yPa$. The commonly used approach based on $G1\u2032\u2264$ 0.95 $GLVR\u2032$ or 0.90 $GLVR\u2032$ underestimates the $\sigma y$ values (Table V). A lack of consistency in these methodologies within the DIW field might unfortunately lead to unreliable metrics based on these values. We defined the flow transition index (dimensionless parameter, $FTI=\sigma f\sigma y$) as the ratio between the stress plateau, or bulk flow stress, divided by the stress value at the onset of the yielding transition.^{27} We found that GO formulations become printable at concentrations that results in FTI $<20$.^{27} All the formulations in this work meet this requirement (Table IV); however, this dimensionless parameter by itself is not enough to guarantee the printability of a formulation, and we advise that other aspects must be considered, such as the magnitude of the storage modulus in the linear regime and the flow stress. This is explained in more detail in Sec. III F.

#### 1. F127 formulations

When active materials (Gr and CNTs powders) are added to F127 stock solution, the behavior shifts from type III (weak strain overshoot) to type I (strain thinning) (Fig. 6). The linear regime becomes considerably narrower for both, with an onset of nonlinearities at $\gamma 0%=0.04$; the storage modulus increases one order of magnitude (up to ∼140 000 Pa); the critical strain, $\gamma c%,$ drops from 11.2 for F127 to 7 for GrCNTsF127; both $\sigma yPa$ and $\sigma fPa$ increase; however, the increase in the yield stress $\sigma yPa$ is more notable [Table IV, Fig. 6(b)]. The extent of the yielding region (for “filled” F127 with active materials) is considerably larger in terms of strain range, but the FTI values are still below 20.

#### 2. GO formulations

When adding active materials (Gr and CNTs powders) to the GO stock solution, two different responses have been observed. *CNTsGO formulation.* The addition of CNTs to GO slightly changes the values of $G1\u2032,\u2009G1\u2033,\u2009\gamma NL%,\gamma c%$, $\sigma y(Pa)$, and $\sigma f(Pa)$, but their behavior is very similar (Table IV, Fig. 7). Both, GO and CNTsGO, exhibit type I behavior,^{61} similar yielding transition, critical strains and plateaus in the stress–strain curve (Fig. 7). *GrGO formulation.* The addition of Gr powders to GO results in a formulation that is not printable (see Sec. III A). The formulation is very stiff, with $G1\u2032$ values of >1 MPa (Table IV), and although it also shows a strain thinning behavior (type I), the critical strain $\gamma c%$ at the crossover point ($G1\u2032,\u2009G1\u2033$) takes place at an even larger strain than GO [∼80%, Table IV, Fig. 7(b)]. The stress–strain curve shows that the stress values for GrGO are considerably higher than for any other formulation [Figs. 7(b) and 7(c)] with a stress value at the start of the yielding transition of $\sigma y(Pa)$ ∼ 463 (Table IV). The extent of the yielding transition also increases due to the onset of nonlinearities at small strains [$\gamma NL%=0.11$], and the increase in the critical strain [$\gamma c%=80$]. Worth noting is the unusual (different to all the other printable samples) stress–strain curve for GrGO [Fig. 7(c)]. It shows a maximum in the meduim-amplitude oscillatory shear (MAOS) region [$\sigma maxPa\u2009\u223c\u20092500;\u2009\gamma max%\u2009\u223c\u20093%]$ at much smaller strains than the critical for this formulation [$\gamma c%=80$]. This extreme behavior can be clearly seen in Fig. 7(c); the stress–strain curve shows a maximum and then drops to a minimum before increasing again at the critical strain. There is not a plateau on stress values that can be correlated with bulk flow. The fact that such a large value of stress is reached at low strain to then drop considerably could be related to slip, shear banding, sample fracture, or microstructure disruption. We cannot confirm or quantify these potential issues, but we can assert that this formulation does not display a smooth yielding transition unlike the other printable samples. Experimental evidence during the printing tests (Sec. III A) suggests that the microstructure becomes disrupted under stress resulting in water segregation at the tip of the nozzle due to filter-press effects; the nozzle gets clogged due to the aggregation of graphite/GO particles that separate from the water.

#### 3. Comparison of Gr formulations

The hybrid formulation containing both additives (F127 and GO) and active materials (Gr and CNTs) shows a very similar behavior to GrF127 (Fig. 8). This network shows a smooth yielding transition unlike the GrGO formulation [Fig. 8(b)]. The hybrid exhibits very similar $\gamma c%\u2009\gamma NL%$, $\sigma y(Pa)$, and $\sigma f(Pa)$ values, and yielding transition compared to those for the GrF127 formulation (Table IV, Fig. 8). These results evidence that the Gr particles and the F127 matrix make the strongest contribution to the rheology of the combined formulation. Although the GO flakes and CNTs in the hybrid formulation do not seem to substantially change the rheological behavior, they do play an important role from a functional perspective.

Plotting LAOS data (from transient or correlation data collection) in different ways [$G1\u2032\u2009Pa$ and $G1\u2033\u2009Pa$ vs $\gamma 0\u2009%$; and $\sigma \u2009(Pa)$ (stress amplitude) vs $\gamma 0$ (strain amplitude) plots] enable us to quantify the values of $\gamma c%$, $\sigma y(Pa)$, and $\sigma f(Pa)$. The stress–strain log–log plot helps identifying the “bulk flow” plateau and any potential issues such as the extreme behavior observed for the GrGO formulation. This valuable and quantitative information might be easily dismissed if only the first harmonic moduli vs strain plot is provided.

### D. Material measures of nonlinearities: $GM\u2032,\u2009GL\u2032,\u2009$*S* (dimensionless index of nonlinearity) and Chebyshev coefficients ($e3e1,\u2009v3v1$)

Using only the first harmonic moduli, it is not possible to quantify and compare the nonlinear behavior of different formulations nor to understand higher harmonics contributions to each of them. This section focuses on the analysis of material measures of nonlinearities for all the formulations in this work. This is to the best of our knowledge, the first time that an FT analysis and quantification of material measures are applied to formulations for DIW. Our analysis focuses on the third harmonic contributions (the most relevant nonlinear contribution) to calculate $e3,\u2009v3,\u2009GM\u2032,GL\u2032,$ and $S$ (equations in Sec. II). First, we quantitatively compare the behavior of formulations prepared with F127 (Fig. 9); we then compare those prepared using GO as the only additive (Fig. 10); and at the end of the section, we provide a comparison of all the formulations that contain graphite prepared either with F127, GO, or a mixture of both (Fig. 11).

#### 1. F127 formulations

F127 stock solution shows an extensive LVR region in the plots for all the material measures of nonlinearities ($G\u2032=GM\u2032=GL,\u2032\u2009S=0,\u2009e3/e1=0,\u2009v3/v1\u22480$, Fig. 9). It then shows an intense transition into the nonlinear regime, as $\gamma 0$ increases into the MAOS and LAOS ($\gamma 0>\gamma c$) regions; the slope on the plot $S$ vs $\gamma 0$ is steep [indicated in Fig. 9(b)]. This suggests that the network of F127 micelles undergoes intense changes in the LAOS region (see L–B plots in Sec. III E). The same trend is observed for the elastic Chebyshev coefficients [determined from the Fourier coefficients, Eq. (2)].^{36} The elastic third harmonic Chebyshev coefficients (relative to the first harmonic) are negative ($e3e1<0$) for all F127 formulations in the nonlinear regime. The trends of the material measures, $e3,\u2009GM\u2032,GL\u2032,$ and $S$, for F127 indicate that the formulations undergo through intracycle strain softening.^{36} Viscous Chebyshev coefficients, $v3v1$ values [Eq. (3)], are *≈*0 in the linear regime; become negative in the SAOS to MAOS transition; show a minimum value in the middle of the MAOS region; and then increase again eventually becoming positive further in the LAOS region [Fig. 9(d)]. This indicates that the viscous component undergoes a transition from a shear thinning behavior ($v3v1<0$),^{36} to a shear thickening behavior ($v3v1>0$)^{36} as $\gamma 0$ increases. The addition of active materials (Gr and CNTs) to F127 stock solution has a considerable impact on nonlinear material measures: $GM\u2032,\u2009GL\u2032$ start to separate from $G1\u2032$ at very small strain values, $\gamma NL%=0.04$, as they enter a MAOS region that spans across a wide range of strain values (Fig. 9). The linear regime (LVR) for these formulations (GrF127 and GrCNTsF127) is dramatically reduced compared to the linear regime for the F127 stock solution (Fig. 9). The active materials change the behavior of the F127 hydrogel, which might no longer be a network of spherical micelles arranged in a cubic lattice.^{46} This re-arrangement of F127 micelles when combined with different solid particles is yet to be studied. All the F127 formulations exhibit $GM\u2032>\u2009GL\u2032$ and $S<0$ in the nonlinear regime, which corresponds to intra-cycle and inter-cycle strain softening.^{35,36} The nonlinear material measure values and trends for GrF127 and CNTsGrF127 formulations almost overlap (Fig. 9); it is clear that Gr powders have the biggest impact in the nonlinear behavior. $S$ and $e3e1$ become negative at very small strains [Figs. 9(b) and 9(c)], with a very gradual as $\gamma 0$ increases through the extensive MAOS and LAOS regions. Linear and nonlinear parameters show that the GrF127 and GrCNTsF127 formulations undergo a more gradual network breakdown into a broader nonlinear regime compared to the F127 formulation base. The viscous Chebyshev coefficients for GrF127 and GrCNTsF127 formulations show a peculiar trend [Fig. 9(d)]; $v3v1$ values gradually evolve from a plateau at slightly negative values ($v3v1\u223c\u22125%$) to positive values in the MAOS region [Fig. 9(d)]. The viscous component of these formulations undergoes a different transition than F127 due to the presence of solid particles; going from a slightly shear thinning behavior to a shear thickening behavior without showing a minimum value on the viscous Chebyshev coefficients [Fig. 9(d)]. These results suggest that the addition of active materials has a different impact on the elastic and viscous contributions to the stress. In this case, the addition of powders to the F127 hydrogel changes the viscous component, while the elastic coefficients follow a similar trend for F127 formulations with and without added particles.

#### 2. GO formulations

The GO stock solution and CNTsGO formulation show almost identical values and trends for all the nonlinear material measures (Fig. 10). However, the not printable GrGO formulation shows clearly differentiated values and trends that are explained in detail in this section. For the printable GO formulations (GO and CNTsGO), the observed trends for nonlinear material measures, $GM\u2032,GL\u2032,\u2009S,\u2009e3,\u2009v3$, are very similar to those described for the F127 formulations containing active materials (discussed in Sec. III D 1). $GM\u2032,\u2009GL\u2032$ separate from $G1\u2032$ in the MAOS region, with $GM\u2032>\u2009GL\u2032$ in the nonlinear regime [Fig. 10(a)]. $S$ and $e3e1$ become negative at very small strains [Figs. 9(b) and 9(c)], with a very gradual as $\gamma 0$ increases through the extensive MAOS and LAOS regions. The not printable GrGO formulations have a narrow LVR, with $G1\u2032,\u2009GM\u2032,$ and $GL\u2032$ values ∼2 orders of magnitude larger (same as for the linear, first harmonic moduli). We found that $G1\u2032,\u2009GM\u2032,\u2009GL\u2032$ start to decrease even within the LVR quantified from the frequency analysis [Fig. 10(a)]. The drop of $G1\u2032,\u2009GM\u2032,$ and $GL\u2032$ values within the LVR for GrGO [Fig. 10(a)] might be associated with viscous nonlinearities [Fig. 10(d)] but not elastic ones, given that $S\u2009\u223c\u20090$ and $e3e1\u223c\u20090$ in the linear regime, while $v3v1<0$. The dimensionless index, $S,$ becomes negative in the MAOS region [Fig. 10(b)], and it reaches a minimum just before the MAOS to LAOS transition, as $\gamma 0\u2192\gamma c$. The elastic Chebyshev coefficient ($e3e1$) follows exactly the same trend [Fig. 10(c)]. The minimum values in both, $S$ and ($e3e1$), suggest that the GrGO formulation undergoes a transition from inter-cycle strain softening to stiffening. This shift has only been observed on this not-printable sample (Fig. 10, and L–B plots and analysis in Sec. III E).

Viscous Chebyshev coefficients for all the GO formulations (printable and non-printable) show the following trend: $v3v1$ values initially fluctuate due to noise at small strains. As $\gamma 0$ increases, they stabilize and take very small negative values ($v3v1<0$, shear thinning) even within the LVR. $v3v1$ values slightly decrease as $\gamma 0$ increases into MAOS, but a change of trend is observed within the MAOS region. The sign of the third harmonic viscous coefficients rapidly becomes positive ($v3v1>0,$ shear thickening),^{36} until they reach a maximum before decreasing again in the LAOS region [Fig. 10(d)]. However similar, the transition observed is more extreme for the not printable GrGO formulation with steeper slopes, with a $v3v1$ minimum value of $\u223c\u221220%$ and a maximum value of $\u223c(80%)$ [Fig. 10(d)]. For this sample, the maximum viscous Chebyshev third harmonic coefficient and the minimum elastic Chebyshev third harmonic coefficient take place at similar strain values $\gamma 0\u223c40%\u201345%$, before reaching the critical strain $\gamma 0=80%$ (Table IV). These results demonstrate that the combination of different additives and active materials changes the elastic and viscous contributions in different and unpredictable ways that could only be revealed through the non-linear analysis. Therefore, future studies will delve into strain decomposition and SPP approaches to improve our understanding of these yield-stress fluids and soft solids.^{13,33,53}

#### 3. Comparison of Gr formulations

A comparison of the material measures of nonlinearities for all the formulations containing graphite powders in this work, including the hybrid formulation, is included in Fig. 11. The trends observed for the hybrid ($S,\u2009e3,$ and $v3$ values) are very similar to all the other printable formulations that contain active materials apart from the not printable GrGO formulation. The hybrid and the GrF127 formulation have very similar linear and non-linear moduli values [Fig. 11(a)], while the not printable graphite formulation (GrGO) is clearly stiffer, with values for these moduli ($G1\u2032,\u2009GM\u2032,$ and $GL\u2032$) an order of magnitude higher.

The quantitative analysis of material measures of nonlinearities for DIW formulations shows that printable formulations with active materials behave in very similar ways, and that the impact of different ingredients (additives and active materials) on nonlinearities can be measured. For example, the printable F127 stock solution shows a “clean” transition from linear to nonlinear behaviors, with an extensive LVR and a rapid breakdown of the microstructure. The other printable formulations in this work (F127 with active materials, GO and CNTsGO) show a very gradual transition into the nonlinear regime unlike F127 by itself. The non-printable GrGO formulation shows some distinctive features in the trends and magnitude of the viscous Chebyshev coefficients. Its behavior shifts from shear thinning to shear thickening as the strain amplitude increases, which suggests that the viscous contributions have an important role in the (not printable) behavior of this sample. Based on the results in this section, and with only one not printable formulation, we cannot yet establish a clear link between nonlinearities and printability. However, we can confidently say that the FT analysis and higher harmonics interpretation provide new insight to better understand yield-stress fluids in DIW, and the role that different additives and active materials play in their rheology.

### E. Lissajous analysis

3D Lissajous spaces [$\sigma t,\u2009\gamma t,\u2009\gamma \u0307(t)$] elastic, viscous, and normalized elastic L–B curves for all the formulations including all the waves are provided in the supplementary material (Figs. S6–S8). Here, we include selected L–B curves in the SAOS, MAOS, and LAOS (F127 and GO formulations in Figs. 13 and 15, respectively) that show differences and similarities between the formulations.

#### 1. F127 formulations

In the *SAOS region* [Fig. 12(a)], the F127 stock solution shows a predominantly elastic response in the LVR up to the onset of nonlinearities ($\gamma NL=2.8%$, Table IV). The widening of elliptical curves in the elastic Lissajous plots evidence that the addition of active materials (Gr and CNTs) to F127 leads to a viscoelastic behavior in the LVR with an increased viscous contribution [Figs. 12(b) and 12(c)]. The ellipse is even wider for the hybrid formulation [Fig. 12(d)]. In the *MAOS region,* the nonlinearities emerge as strain increases and non-elliptical distortions appear in the L–B plot. The large and minimum storage moduli $GM\u2032$ and $GL\u2032$ can, in principle, be determined directly from the L–B plots [Fig. 12(a), MAOS].

From a careful look on the L–B shapes, it strikes that $GL\u2032>GM\u2032$ for all the formulations in the MAOS region, which contradicts our Chebyshev analysis (in Sec. III D). To confirm this observation, the elastic nonlinear moduli for the F127 stock solution have been calculated and compared using the two methods: (1) from the Chebyshev analysis including higher harmonics up to n = 11, and (2) from the local slopes in the L–B curves at zero strain (tangent modulus, $GM\u2032$), and at maximum strain (secant modulus, $GL\u2032$). The results of this analysis confirm the discrepancy (Fig. 13). The values calculated using the two methods match within the LVR; however, they divert in the nonlinear regime, thus confirming that $GL\u2032>GM\u2032$ based on the L–B plots, which is associated with strain stiffening behavior within a LAOS cycle,^{36} while the Chebyshev analysis including harmonics up to *n* = 11 confirms that $e3$ takes negative values [Fig. 13(a)] and that $GL\u2032<GM\u2032$ for all the strains in the sweep [Figs. 13(c) and 13(d)], which are associated with strain softening. This rather confusing contradiction or “paradox” has been previously discussed in the literature, and that this can be explained by the sensitivity of local measures of elasticity to the rate of the deformations, which is not constant within a LAOS cycle.^{38,57} The fact is that all the formulations are clearly strain softening as the strain amplitude increases in the sweep, and that $G1\u2032,\u2009GM\u2032,$ and $GL\u2032$ all considerably decrease in LAOS, up to ∼3 or more orders of magnitude.

#### 2. LAOS region

The L–B curves for F127 [Fig. 12(a), LAOS] get close to a perfect square (when axes are normalized, Fig. S7 in the supplementary material) in the nonlinear regime. Square shapes in elastic Lissajous curves are associated with plasticity,^{38} and a perfect plastic response would correspond to the maximum possible dissipated energy for a given strain amplitude $\gamma 0$ and maximum stress $\sigma max$. In the L–B curve for F127 at the largest strain amplitude in the sweep [Fig. 12(a), LAOS], we can also point out the yield point that shows the transition from purely elastic behavior, to yield into viscoelastic and viscous flow at the top of the curve (plateau at the maximum stress $\sigma max$). The perfect plastic dissipation ratio ($\varphi $) defined in Sec. II [Eq. (7)] is said to be a well-defined material measure with almost no sensitivity to the data acquisition rate (or rate of deformations).^{38} Here, we use this metric to compare the formulations in this work and to quantify how close or far they are from the perfect plastic behavior. F127 shows an almost purely elastic behavior in SAOS ($\varphi \u2009\u223c\u20090$, Fig. 15) and gets close to an ideal plastic behavior in LAOS as $\varphi \u21921$. Xanthan gum and drilling fluids are examples of yield-stress fluids that exhibit some similarities with the F127 stock solution.^{38} However, when active materials are added [Figs. 12(b) and 12(c)], the curves in the LAOS region change shape, and the transition between elastic deformation and viscoelastic flow is not easy to distinguish within a LAOS cycle, because the L–B shape becomes smoother with softer edges in the LAOS region [Figs. 12(b) and 12(d), LAOS]. The dissipation ratio trends for F127 formulations confirm that the addition of active materials increases the viscous contributions in the SAOS region [$\varphi \u2009\u223c\u20090.3$, Fig. 15(a)]. The dissipation ratios in LAOS for GrF127, CNTsGrF127, and hybrid formulations show values in the range between $\varphi \u2009\u223c$ 0.7 and 0.8 [Fig. 15(a)], which is associated with Newtonian behavior instead of ideal plastic. These results suggest that the pure F127 stock solution can be classified as an elastoplastic material, while “filled” F127 formulations with carbon materials are viscoelastic.

#### 3. GO formulations

The L–B plots for the printable GO formulations [Figs. 14(a) and 14(c)] show similar shapes and trends than the filled F127 formulations with active materials (for all three SAOS, MAOS, and LAOS regimes), but the L–B plots for the non-printable formulation, GrGO [Fig. 14(b)], has very distinctive shapes in the LAOS region.

In the *SAOS region,* GO stock solution and CNTsGO formulation show a predominantly elastic behavior based on both the L–B plots [Figs. 14(a) and 14(c)] and the values of dissipation ratio [$\varphi <0.1$, Fig. 15(b)]. The L–B plots show elliptical distortions at very small strain values for the non-printable GrGO formulation [Fig. 14(b)]. In the *MAOS* and *LAOS regions*, the differences between printable and non-printable GO formulations become evident. The GrGO formulation shows extreme distortions [Fig. 14(b)], with upturns at the maximum strain within a LAOS cycle. It is also possible to identify the peak of the stress amplitude ($\sigma max)$, which has been discussed in detail in Sec. III C. The values of the dissipation ratio, $\varphi $, for this formulation show a non-monotonic trend unlike the printable formulations [Fig. 15(b)]. In SAOS, $\varphi $ varies arbitrary between 0 and 0.2; it then increases in MAOS reaching a first maximum value $\varphi \u2009\u223c\u20090.4$ at 5% strain amplitude ($\gamma 0$); it then goes through a minimum and increases again up to $\varphi \u2009\u223c\u20090.5$ at $\gamma 0$ = 100%.

The visual representations of LAOS results in L–B curves provide a very useful tool to identify nonlinear behaviors of printable (or not) formulations. These plots can be easily generated from a standard LAOS test without needing to delve into a quantitative nonlinear analysis. DIW users can benefit from monitoring and reporting these curves, to help with the analysis and discussion of the yielding transition, and to compare formulations with different compositions. Even in the absence of frequency analysis or quantitative measures of higher harmonics contributions, Lissajous plots will help researchers to better understand their formulations and identify potential issues in their microstructure and printability. The determination of the dissipation ratio combined with the Lissajous analysis complements the findings from the frequency analysis (Sec. III B), the stress–strain plots (Sec. III C), and the quantitative materials measures of nonlinearities (Sec. III D). This “well-behaved” scalar measure^{38} with almost no sensitivity to the acquisition rate seems capable of differentiating printable and non-printable formulations. Although a clear link between the dissipation ratio and printability cannot be drawn quite yet given the reduced number of formulations studied here, $\varphi $ is a promising metric to characterize the behavior of DIW formulations.

### F. Printability maps using Ashby-type material measures charts

In Secs. III B–III E, we provide the quantitative analysis and qualitative visual representations of LAOS results for a library of soft solids designed for 3D printing of conductive 3D architectures. We have described the concept of printability and the printability of our formulations in Secs. I and III A, although an automatized system to objectively judge “print quality” has not been used, and therefore, this judgment might be subjective. However, in this relatively broad definition we also must account for the wide range of applications involved and that printability requirements and maps need to be tuned accordingly. This determines the level of control and accuracy to realize complex features in 3D shapes. Bioprinting applications might require moderate $G\u2032,$ and low $\sigma f$ values to guarantee cell survival and proliferation,^{8–10} while other applications might require an accurate control of fine features, such as unsupported overhangs, which need higher $G\u2032,$ and $\sigma f$ values. The map here discussed focuses on formulations of carbon-based and other 2D materials for energy applications. We have previously described the printability of soft materials based on the following bulk rheological parameters: $G\u2032,\u2009\sigma f,$ and $FTI$.^{27} The storage modulus ($G\u2032$) and the yield (or flow) stress ($\sigma y*$) are frequently used to describe printability.^{31,34,62,63} However, note the asterisk to flag the lack of consistency within the DIW field on the yield (or flow) stress determination. For example, $\sigma y*$ could be determined from fitting continuous shear (flow ramps) data; from the stress determination at the end of the LVR; or the stress at the solid-to-liquid transition ($G\u2032=G\u2033$). This leads to inconsistencies and lack of consensus when discussing printability. Similarly to Ashby's materials selection charts,^{39} we can create maps that compile linear material measures $G\u2032$, vs $\sigma y*$ on the coordinate axes for different formulations, that can be obtained from a standard LAOS test without the need of doing an FT analysis. We provide here Ashby-type printability maps using the results from this work [Fig. 16(a)], and for context, we populate this map with data from other formulations found in published literature [Fig. 16(b)]. A clear correlation between nonlinearities and printability cannot be concluded from the FT analysis performed on a relatively small set of samples (only just one of them being not printable). What we can confidently conclude from our analysis though, is that the trends of nonlinear parameters for all the printable formulations are the same. Our long-term aim is to expand our analysis to a larger set of samples in DIW for different applications and expand these printability maps to include linear and nonlinear material measures, and scalar metrics such as the dissipation ratio. We will also delve into additional tools, such as the sequence of physical processes (SSP) developed by Rogers *et al.*,^{50–52} and however, this is currently beyond the scope of this manuscript.

The printability map, based on linear material measures, for the formulations studied in this work shows that all the printable formulations are contained within a small region [highlighted in turquoise, Fig. 16(a)], while the non-printable GrGO formulation is far from it, due to its high stiffness and the stress peak on the stress–strain curve (Sec. III C). The map also includes dashed lines delimiting the region of interest in which our formulations fall: two vertical lines for flow stress limits ($\sigma f=250\u2009Pa$ and $\sigma f=1500\u2009Pa$) chosen based on our empirical evidence, and a diagonal line that represents the figure of merit ($FoM=G\u2032\sigma f=20$).^{62} The printability region is not “rigid,” and its boundaries can be blurred: softer formulations ($G\u2032\u2193,\u2009\sigma f\u2193$) located toward the left on the map and close to the diagonal (area highlighted in light color) might be suitable for 3D printing of parts with limited complexity. As we move toward the right and upward from the dashed diagonal on the map, we find stiffer formulations ($G\u2032\u2191,\u2009\sigma f\u2191$, area highlighted in darker color) that provide better print quality.

For example, the F127 stock solution can be 3D printed to make relatively simple parts with vertical walls [images in Fig. 16(a)], while the hybrid formulation can also be used to build complex parts with fine features [images in Fig. 16(a)]. It must also be considered that the density and concentration of the materials in the formulation, and whether the printing takes place in air, inside another gel, or an oil bath, will affect the required $G\u2032,\u2009\sigma f$ values for that formulation to be printable.^{64} We find some limitations to build reliable maps in context with data from the literature, because the quality of the data collected is often uncertain, or the analysis is not explicit [Fig. 16(b)]. Despite the barriers, we find that most of the formulations labeled as printable in the literature here included fall within the printable regions discussed above. These maps can be a very useful tool in formulation design if they were carefully populated with reliable data. However, we affirm that they should not be used in isolation, because this could lead to an oversimplification of the analysis and overlooking signs of unusual behaviors.

## IV. CONCLUSIONS

The rheology of yield-stress fluids and studies of their LAOS behaviors are active and long-standing research fields, while DIW or 3D printing of soft matter is an emerging research area that is quickly expanding. Here, we aimed to build a bridge between more specialized rheology studies on yield-stress fluids, with the growing diverse community exploiting DIW in different applications. This work provides to the multidisciplinary community in DIW, guidance to improve data collection, and tools to perform deeper analysis of LAOS tests (including visual and qualitative analysis of Lissajous–Bowditch curves; yielding transition quantification in stress–strain curves; frequency analysis and quantification of higher harmonics, as well as material measures of nonlinearities). This work has focused on LAOS of soft solids designed for 3D printing carbon-based conductive structures. Our library is a small but fair representation of the different materials and formulation approaches currently being developed, and unintentionally has provided us with an opportunity to compare printable and non-printable formulations. From the results organized in different sections, we can draw the following conclusions summarizing our findings.

FT rheology of raw data signals enables us to quantify nonlinearities and to create 3D maps to visualize the onset and rising of higher harmonics contributions. We have shown how 2D projections of these maps can be used to carry out quantitative comparisons of different formulations and to identify the role that different additives and active materials play in the yielding transition. Quantitative material measures of nonlinearities (Fourier coefficients ($G1\u2032,\u2009G1\u2033,G3\u2032$, $G3\u2033,$ etc.); Chebyshev coefficients ($e1,\u2009v1,e3$, $e3,$ etc.); elastic moduli $(GM\u2032,\u2009GL\u2032);$ and dimensionless index of nonlinearity ($S$) have enabled us to explore the higher harmonics contributions and how they change for different additives and active materials. Although all the formulations in this work have shown broadly similar trends (most of them showed type I behavior), each of them has characteristic values for material measures in the linear and nonlinear regime. We have objectively determined: the onset of nonlinear behaviors (end of LVR, $\gamma NL$) from harmonics maps; and the yielding transition and bulk flow stress ($\sigma f$) from stress–strain graphs. Finding the limits of the LVR for DIW formulations can be extremely subjective, and its end is sometimes delimited considering a drop to ∼90% of the at rest storage modulus, $GLVR\u2032$. The FT analysis enables to precisely quantify the onset of nonlinear behaviors and to evaluate how this is affected by the addition of different active materials and additives. This analysis will be very useful to those working in bespoke additives and formulation design, for example, in the synthesis of new generations of responsive surfactants.

We have encountered contradicting results between the Chebyshev analysis and the determination of $GM\u2032,\u2009GL\u2032$ from local slopes in the L–B curves that can be explained due to the sensitivity of local measures of elasticity to the data acquisition rate. An alternative scalar metric, the dissipation ratio, $\varphi $ (which is less sensitive to the rate of deformations), is a promising parameter to characterize the behavior of complex fluids designed for DIW.

The stress–strain log –log plots and the Lissajous–Bowditch curves displayed singular features for the non-printable formulation (GrGO) that clearly differentiate it from printable ones. This formulation reached an extremely high stress peak or maximum within the yielding transition instead of the plateau associated with bulk flow that printable formulations displayed. This undesired behavior is clearly seen in the stress–strain plot, in the non-monotonic trend of the perfect plastic dissipation ratio, $\varphi $, and the L–B plots in the LAOS region. Our analysis and results flagged issues that would have been overlooked if the analysis was based only on the first harmonic moduli ($G\u2032,\u2009G\u2033$ vs $\gamma 0$) plot. The viscous Chebyshev coefficients have also shown an extreme behavior for the non-printable formulation. We have demonstrated the value of performing a deeper analysis of the raw data signals, quantitatively (through frequency analysis and quantifying material measures of higher harmonics contributions) and qualitatively (through Lissajous analysis) to achieve a deeper understanding of printable and non-printable formulations.

Ashby-type charts based on linear material measures, $G\u2032$ vs $\sigma f$ on the coordinate axes, are very useful printability maps to guide formulation design and pave the wave for common protocols and standards in DIW. In this map, focused on formulations for energy applications, we have delimited our printable region using the results from this work. We have expanded its boundaries using data from the literature and discussed the needs of different applications. Printable formulations can range from “soft” for the printing of simpler 3D shapes with limited complexity (for example in bioprinting applications) to “stiff” for the accurate printing of complex shapes with fine features. We would like to emphasize the importance of using these maps alongside other quantitative measures (stress–strain log–log plots) and/or visual representations (Lissajous analysis) of LAOS results to avoid overlooking the signs that we have associated with poor printing performance. Although the link between nonlinear materials measures and printability does not clearly emerge yet from the analysis of a small set of formulations, we have demonstrated that the combination of quantitative material measures, visual Lissajous–Bowditch representations, and Ashby-type charts are a powerful approach to fully characterize complex fluids in 3D printing, and to guide the design of new additives and formulations.

## SUPPLEMENTARY MATERIAL

See the supplementary material for additional information in data collection, additional figures, and discussion to support Secs. III B (stress waves and FT 3D maps for all the samples are provided) and III E (including 3D Lissajous spaces, elastic, and viscous projections with all the waves for each of the formulations).

## ACKNOWLEDGMENTS

This work was supported by a UKRI Future Leaders Fellowship MR/V021117/1.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Esther Garcia-Tunon:** Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). **Rishav Agrawal:** Formal analysis (supporting); Validation (supporting); Writing – review & editing (supporting). **Bin Ling:** Methodology (supporting). **David JC Dennis:** Data curation (supporting); Formal analysis (supporting); Software (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

_{3}C

_{2}T

_{x}(MXene) ink for fabrication of micro-supercapacitors with ultra-high energy densities