Direct Writing Unclonable Watermarks with an Electrochemical Jet

Counterfeit parts result in significant losses per annum and are often dangerous, therefore they represent a serious concern for manufacturers and end users alike. Easily written but unclonable watermarks undermine the proposition of the counterfeiter. Here, a rapid electrochemical jet engraving routine is presented to encode robust materials with self‐organized dendritic structures at length scales that can be imaged with a smartphone. Surface defects act as stochastically distributed seeds from which discrete pitting events can be propagated by translating the electrochemical field. While the vascular pathways can be directly written at the macro scale, the formation and propagation of microscale dendritic arms is chaotic, caused by the implicit randomness of the defect seeds and the supply of ions to the surface. The latter is confounded by random perturbations in the flow condition. Each engraved dendrite is unique, stable at high temperature (>500 °C) and can be subjected to rapid image recognition to allow individual mark identification at any point during part production and delivery, or through part lifetime.


Introduction
Complex self-organized flow architectures are ubiquitous in natural systems across length scales (Figure 1a-e). In biological systems, evolutionary pressure has guided the generation of optimized mass transport networks through natural selection. Geological features like river basins [1] and even celestial bodies are postulated to undergo similar morphodynamic evolution, [2] with constructal theory stating: "for a finite flow system to persist in time, it must evolve in such a way to give easier access to the imposed (global) currents that flow through it". [3,4] Our attempts to mimic natural self-organized flow structures have physical unclonable functions (PUFs). While discrete to watermarks, PUFs have features with inherent randomness created by non-deterministic or pseudo random processes. [18][19][20] PUFs with implicit randomness (arising from variations during manufacture) cannot be replicated and provide useful entropy sources generating millions of unique keys. This has applications wherever a cryptographic challenge-response is required. While generally associated with electronic systems, PUFs can be incorporated into consumer products, documents, and even pharmaceuticals to deter the counterfeiter. [18,[21][22][23][24] PUFs exploiting the unique arrangements particles make during drying or coating have been used to authenticate surfaces. Bae et al. created polymeric particles with random wrinkles within a silica film leading to fingerprint-like minutiae imageable with a smartphone camera. [25] Liu et al. to inkjet printed CdSe quantum dots to direct write PUF tags, [26] fabricating features visible to the naked eye upon UV irradiation. Carro-Temboury et al. exploited the discrete luminescence behaviour of three lanthanide ions (Eu 3+ , Tb 3+ , Dy 3+ ) to create unique patterns upon irradiation at different visible light wavelengths. [27] Gu et al.   gap-enhanced Raman tags by drop-casting nanoparticle dispersions, reporting extremely high encoding capacities (e.g., > 10 15 051 ) when applying multiple nanoparticle types, requiring Raman microscopic analysis. [28] While these represent important advances to the state-of-theart, several issues must be addressed to apply stochastic watermarks to engineering parts. Particle tags and films may not be compatible with high thermomechanical loads or hard-wearing applications. A requirement for microscopic analysis or adaption of conventional imaging systems complicates end-user authentication and reduces practicality. The bespoke nature of certain taggants, like nanoparticles with embedded Raman active material, challenges scale-up and deployment as a practical anti-counterfeiting routine. Finally, surfaces decorated with nanoparticles, which may penetrate the dermal barrier, are unlikely to be adopted in human-facing applications such as consumer goods, particularly those bearing heavy metals like Cd. There remains a requirement for a codification system that is robust, easy to write (and scalable), easy to read, suitable for human interaction, while being impossible to duplicate.
Here, we introduce a method through which unique watermarks can be directly written onto engineering materials as a low-cost, ambient condition jet electrochemical etching routine [29,30] (Figure 1f). Electrochemical marking can be applied to metals and alloys, as well as semiconductors. These classes dominate material selection in high-value consumer goods and engineering applications, particularly where robustness is a requirement. Within electrochemical systems, mass flux occurs by diffusion, electromigration, and convection (advection) mechanisms ( Figure 1g). These are the 'flow' mechanisms that can be exploited to self-organize vascular networks when the field is translated over a surface (Figure 1h,i) at mm-length scales (following the deterministic toolpath) visible to the naked eye that can be imaged with a smartphone (Figure 1j). The propagation of these networks is directed partially by turbulent flow and diffusion and initiated at stochastically distributed surface defect sites (pseudo-random seeds); there is implicit randomness in fabrication. Thus, unclonable watermarks can be applied to robust engineering materials using a factory floor compatible machine tool. While high-value parts are one target application of this phenomenon, the findings from this study can be also be applied to consumer goods, luxury products, or in any application where there is a desire for authentication and traceability.

Engraving Dendrites Beyond the Diffusion Limit
The anisotropic growth of electrochemical corrosion pits leads to radially directed surface features as interaction time is increased (Figure 2). Initially, relatively small interaction times (0.25 s, Figure 2ai) lead to the formation of discrete isotropic corrosion pits in the interaction region beneath the nozzle (current density, J = 100 A cm −2 ). This is distinct to hemispherical erosion and ablation observed in other energy beam processes [31][32][33] including lasers and electron beams, as well as electrolyte jets in other electrolyte-material systems. [34] The distribution of these initial pits appears stochastic within the interaction region. Their distribution, however, is not truly random, but indicates that initial breakdown is likely to depend on the local integrity and availability of defects in the native oxide layer. It is these initial pits that act as the seeds for subsequent propagation. As the interaction time is increased to 0.50 and 1.00 s (Figure 2aii,iii), morphological anisotropy increases, leading to the development of directional etching channels influenced by the radial flow of the wall jet region and the current density distribution. [35] The isotropic to anisotropic dependence on timescale is consistent with multi-mechanism flow phenomena. [36] Perceptually, this is analogous to self-organization phenomena like crystallization. In heterogeneous nucleation, particles initially aggregate around foreign matter, before growth continues along the crystal habit or opposite to the heat flow, forming the familiar symmetries associated with snowflakes and dendrites (Figure 1d,e).
As charge is transferred, pitting likely initiates at surface regions associated with local imperfections in the oxide layer, [37] like scratches (Figure 2b), at grain boundaries, or where the surface is under local residual stress. [38] Ti is likely to oxidize to Ti(IV) under these conditions [30] according to Equation (1), the rate of which is sensitive to the surface concentration of iodide, [I¯]. Initially, this is the bulk iodide concentration, and pits can develop rapidly. In addition to metal liberation, iodide is easily oxidized at the anode, [39] electrodepositing molecular iodine according to Equation (2), depleting [I¯] and reducing the rate of Equation (1). Charge is also transferred through anodization, simplified in Equation (3), dehydrating the iodine-rich film adjacent to the anode. As these reactions occur, a poorly soluble iodine film forms over the anode, the surface concentration of [I¯] tends toward zero, wherein the limiting current condition is exceeded, and diffusion becomes the dominant mechanism through which iodide is supplied to the surface to remove material. This film is thinned by forced convection from the impinging jet and diffusion, in addition to comproportionation to form triiodide [40] (Equation 4), driven by the high solubility of triiodide in water. The film reaches an equilibrium thickness determined by the local fluid velocity [39] and the relative rates of these reactions.
Material removal in this equilibrium system is partially reliant on the diffusion of iodide through the film to the anode, governed by Brownian motion. Assuming a passivating oxide layer acting as a partial mask, iodide cannot partake in Equation (1) except at a pre-existing pit or oxide layer defect. Iodide performs a "random walk" through the film and the anode until it reaches a pre-existing pit. This leads to pit morphologies typical of diffusion-limited aggregation (DLA) (Figure 2c), consistent with other diffusion-limited electrochemical phenomena, like electrodeposition. [41,42] It is notable that material removal is inhibited in the central region at machining times of 1.00 s, despite the center being the point of greatest current density. [35] This can be explained by screening of the central region toward diffusing iodide by the propagating branches. [43,44]  + + . This supports the removal reaction being limited during a given exposure, for example through a diffusion-limited mechanism (Figure 2c). At short exposure times (<0.50 s) removal rates are greater than when the exposure time is lengthened. It is assumed that the iodide supply to defect sites or pre-existing pits dictates the Ti electrolysis rate. Where surface iodide concentration reduces, the rate of Equation (1) reduces along with the charge-specific removal rate. Increasing applied current density also reduces the charge-specific removal rate, which is intuitive given the rate of iodide depletion will be increased.
Shape descriptors like circularity describe how morphologies evolve with exposure time at different applied currents. Circularity ( Figure S1, Supporting Information) reduces as exposure time increases as the pits radially propagate (Figure 2e). Circularity is retained more effectively at the lowest current density (50 A cm −2 ), compared with the higher current densities tested. This results from two factors, i) the reduction in iodide depletion rate at low current densities, and ii) from fewer breakdown sites resulting from each given interaction (Figure 2f). The latter is a  1.00 s (J = 100 A cm −2 , νj = 23.4 m s −1 ). As exposure time increases, pit morphologies become more anisotropic and radially directed from the center of jet impingement. b) Initially, the anode surface concentration of iodide is high, initial corrosion pitting locations likely depends on local surface oxide integrity, where pitting propensity is greater at defects like scratches shown in the electron micrograph. Iodide oxidation also occurs leading to the formation of an anodic iodine film, depleting the surface concentration of iodide. c) As iodide surface concentration tends to zero, a diffusion limited steady state is achieved, and material removal is reliant on reactant species diffusion to the surface, leading to dendritic structures growing radially from the central impingement region. d) Charge specific removal decreases with interaction time, after that an apparent steady state is achieved. e) Circularity of pits decreases with exposure time as radial directionality increases. f) Pitting frequency decreases with exposure time and applied current density. g) Pits can be propagated by translating the jet (J = 75 A cm −2 , νf = 0. function of the reduced charge quantity and the lower potential applied to the anode to fabricate the marks for given exposure times at 50 A cm −2 . Fewer pitting sites within the interaction region will reduce screening and make combinations of discrete pits rarer. Translating the electrolyte jet imposes a directional bias to pit propagation, retaining the directional morphologies (Figure 2g), where bifurcating and interconnected dendritic structures are generated along the length of the toolpath (J = 75 A cm −2 , feed rate νf = 0.3 mm s −1 , jet velocity νj = 24.3 m s −1 ). These patterns are akin to fluid displacement patterns. [45,46] The depths of the branches appear largely independent of the local current density distribution, which is consistent with the removal mechanism being controlled by diffusion of iodide, where attack is more likely at the edge of a given branch compared with the basal area. Finer and shallower features in the vicinity of the translation finish point result from the reduced interaction time and local charge passed. While removal is defined by the deterministic translation of the jet, channel architectures are chaotic.
The micrograph ( Figure 2h) shows a wall section of Ti between two micro-branches. The tilted view (65°) indicates undercutting of the surface oxide, also shown in Figure 2i, which reveals an etched microscale surface characterized by protruding β-phase (white). This is typical of electrochemical corrosion in dual-phase Ti-alloys. [47] Undercutting has previously been reported when electrochemically processing Tialloys, where the TiO 2 film acts as a mask. [47,48] This is consistent with the prevalence of side wall attack, which would be expected from a mechanism reliant on diffusion and migration of the electroactive species across the anode surface. Further evidence to support the theory that propagation is predominantly controlled by the growth of pre-existing features, rather than the nucleation of new corrosion pits ahead of the channel, is shown in Figure 2j. Here, the leading edge of the channel appears with no observable adjacent pitting. In this mechanism, the native oxide acts as a mask through which charge transfer is focused once the initial pits are generated. [30] Potential traces acquired during exposure allow pit evolution to be understood. Figure 3a,b show binary topographic masks fabricated at two different current densities (50 and 150 A cm −2 ) for equivalent features marked with the film-forming iodide electrolyte compared with another chaotropic halide anion of similar electrolytic conductivity (NaCl, 147.7 ± 0.1 mS cm −1 at 21.8 °C). Chloride electrolytes used to electrochemically machine Ti-alloys, [47] allowing high material removal efficiencies at high current densities (>150 A cm −2 ), [30] indicated by the large pits compared with iodide. Potential traces acquired when marking these features are shown for both NaI (black series) and NaCl (red series), where behaviour depends on electrolyte, applied current, and exposure time. At low current density (50 A cm −2 ), the potential falls during each interaction with the iodide electrolyte, which is most noticeable for the longer exposure times (2.00 and 5.00 s). This is likely caused by the high charge transfer resistance over the passive TiO 2 film. During each interaction, pits propagate laterally as material is removed and the contribution from charge transfer resistance drops.
Contrarily, with chloride, potential increases during each interaction resulting from an increase in the ohmic drop as material is removed. [49] This is noticeable at high current density (150 A cm −2 ) where more Ti is removed. With iodide, a fraction of the charge passed is consumed by iodide oxidation rather than Ti electrolysis and the ohmic drop increase smaller. This is shown in the reconstructions of typical marks in Figure 3c, along with extracted profiles for both electrolytes (5.00 exposure times). Marks fabricated with chloride have larger areas at both current densities and are much deeper at high current density. At high current density, the potential decrease during interaction with iodide is less pronounced, correlated to increased breakdown and greater pitted area at higher currents. Thus, the contribution of this exposed area to the charge transfer resistance over the surface is smaller. Pitted area is inversely correlated to the final potential (measured from the last 150 ms, Figure 3d) and is most acute at low current density (50 A cm −2 ), where the pitted area is generally small.

Flow Directed, but Unique?
Evidencing uniqueness presents an interesting philosophical and logical challenge. Proving that sequences of independent and interdependent stochastic events will always lead to a distinct result cannot be guaranteed, but the likelihood can be shown to be extremely low. To address this, the sequence of events and associated physical phenomena resulting in a pattern must be considered, where pattern formation occurs in three general steps: i) initiation, ii) propagation, and iii) termination.
There is a high probability that an initiation point will occur within two nozzle diameters (2d N ) of the impingement zone (Figure 4), but the location cannot be predicted. Consider the polycrystalline substrate with a nonuniform native oxide layer. While initial breakdown can occur at identifiable features like surface scratches (Figure 2b), more often this does not appear to be the case. Breakdown also occurs at stochastically distributed oxide defects and vacancies [50] and pitting is also influenced by the presence and character of grain boundaries. Taken together, these provide a rich entropy source at smaller orders of magnitude than the deterministic length scale. Typical surface grain boundary densities between 3-5 mm mm −2 would be expected from this substrate in the rolled supply condition (grain size 2-10 µm, Figure 2i), while oxide defect densities between 10 17 -10 19 mm 2 might be expected depending on the condition of the oxide film. [50,51] Considering most initial pits form within 2d N (0.82 mm 2 ) of impingement, predicting the precise initiation point of 10 s of pits (Figure 2f) becomes an unrealistic proposition requiring an intimate knowledge of the surface that far exceeds current capability.
Propagation depends on the flow of reacting ions to the surface and a pre-existing pit. Flow mechanisms like diffusion and forced convection within the turbulent flow regime cannot be predicted to the extent that local concentrations can be accurately deduced over a surface. One can consider propagation as an array like that observed in a Monte Carlo simulation where the precise location of material removal is dependent on preceding removal, but where future events cannot be reliably repeated. Importantly, the formation of localized limits terminates propagation, which is also commonly observed in other self-organization phenomena like crystallization. As such, no two interactions should result in features with identical morphologies. The challenge then orients toward deploying these watermarks in a manner that is practically deployable for verification purposes (they should be easy to write and easy to read).
The high optical contrast of the pitted regions relative to the top surface makes the marks particularly amenable to basic image processing operations such as thresholding and binarization that are not computationally expensive. Individual marks were processed at identical jetting parameters (100 A cm −2 , 0.75 s dwell time) to demonstrate individuality. Figure 4a shows a sample piece (20 × 20 mm) with multiple marks with an example image processing flow from a raw photograph of a given mark to feature extraction. Features were extracted from raw images following a 'skeletonization' morphological filter, [52] from which minutiae such as relative branch and termination point positions can be determined. This routine was inspired by conventional fingerprint identification approaches. Figure 4b shows the radial distances, R, of branch and termination points from the respective image centroids aggregated from 60 discrete marks, where branch points are clustered closer to the center of the interaction region (median, 0.20 mm) than the termination points (median, 0.36 mm). Intuitively, termination points are more numerous (n = 3993) than branch    Figure 4j were generated from the polar and radial data strings extracted from 60 individual patterns. In each histogram, the Weibull distribution of the data is presented. points (n = 527), where the branching ratio will depend on the jetting parameters, and every branching instance will result in at least n + 1 termination points. For these jetting parameters, most identified branches fall within two nozzle radii (95th percentile, 0.47 mm), while 15% of the termination points fall outside two nozzle radii (85th percentile, 0.50 mm). 95% of the branch and termination points fall within 1.14 mm 2 surface area processed in 0.75 s. Figure 4c shows binary masks extracted from photographs of these 60 marks. In each case, the specific mark is unique. Pitted areas for 40 marks processed on the same sample piece (Figure 4d) are clustered between 0.113-0.151 mm 2 (5 -95th percentile), indicating areal similarity despite morphological variance. Figure 4e shows the normalized superimposition of the 60 marks within rows 1-6. Individual marks possess a stochastic pattern defined by the substrate and the process; this is shown clearly in mark 1-vi imposed over two approximately parallel surface scratches. However, aggregated material removal follows a distribution correlating with the applied current density over the surface. [29,35] Extracted features can be applied to identify individual marks using image recognition processes. [26] In this study, we present a rapid approach to acquire data strings from engraved marks using relative locations of extracted termination points for comparison through computationally inexpensive logical operations. Each termination point within a radial mark can be described by a radial distance, R, and a polar angle, θ, from the image centroid, shown in Figure 4f from the specific mark shown in Figure 4a. The radial distance (in pixels) and the polar angle (degrees) of each termination point can be written to a bit string, where a positive bit corresponds to one or more termination points at a given sampling interval (pixels or degrees) shown in Figure 4g,h (extracted from Figure 4f). This approach is robust to image rotation, nevertheless alignment can be assured in a trivial manner simply by engraving three adjacent marks (consider QR code finder patterns), requiring a small cumulative exposure time (2.25 s) at these marking parameters. As indicated in Figure 4b, the radial data string is clustered away from the image center (0-pixel point) and at the image extremes. The polar string is more evenly distributed, which is intuitive considering there is no apparent bias in propagation direction over the 60 marks (Figure 4e). In the polar string, positive bits cluster together as a given pit propagates radially.
n Terminations <n bits within a given bit string, i.e., a positive bit corresponding to a termination point is a low probability (highinformation) event, shown in the entropy plots (Figure 4i) for the radial and polar strings respectively. High information and low entropy within the data strings is grounded in the observation that there are relatively few termination points, and the marks are localized within an area determined by the jet interaction. This betrays the fact the marks have both deterministic and random factors. It should also be noted that the measurement of entropy is sensitive to the sampling interval [53] (a 1:1 bit to pixel or polar angle ratio was applied here). Data compression algorithms have been applied to address similar issues with low bit densities in PUFs for example. [24] This also applies to the Hamming distance, defined as the number of bits that require inversion, or 'flipping', to make one string identical to another. Figure 4j shows the distribution of Hamming distances for the 60 radial and polar data strings. Given the polar data strings have a rotational dependence, the Hamming distance between polar strings was calculated by iteratively shifting one string relative to another in a circular manner through each sampling interval (i.e., a rotation), and where data reported is the minimum Hamming distance acquired from this operation. The mean Hamming distance for the 60 radial strings is 109 bits (90 -129 bits 5 -95th percentile), and the mean of the minimum polar Hamming distances is 78 bits (64 -90 bits 5 -95th percentile). An additional encoding step considering branch point locations can be adopted in a trivial manner, however the theoretical encoding capacity considering solely the two data strings associated with the termination points indicates that one mark can be readily discriminated from another.
The similarity between extracted radial and polar data strings acquired from the 60 discrete marks is shown in the correlation matrices in Figure 4k,l, respectively. Aggregated data for the radial and polar data strings are given in Figures S2 and S3, (Supporting Information) respectively. The mean of the nonidentical values is 0.043 ± 0.042 (min = −0.071, max = 0.211) for the radial strings and 0.002 ± 0.055 (min = −0.167, max = 0.169) for the polar strings, indicating that marks are poorly correlated with one another. The marginally more positive correlation between the radial data strings is likely related to their distribution around the image centroids. Logical cross correlation between all 812-bit radial strings took 0.555 s, and all 360-bit polar strings took 0.519 s using a standard desktop PC (Intel Core i3-7100 @ 3.90 GHz, 8 GB RAM) running the inbuilt 'corr' linear correlation coefficient MATLAB function.

Controlling Morphodynamic Evolution
Stochastic pit morphologies are retained during jet translation, allowing unclonable patterns to be directly written onto a surface. Figure 5a-c shows reconstructions of typical surfaces marked between 75-150 A cm −2 . Electrolyte flow velocity (ν j ), translation speed (feed rate, ν f ) and applied current density (J) affect the stability of pits as they propagate across a surface. These are the key characteristics, along with the nozzle diameter, that dictate the scale of features in conventional electrochemical jet processing. [54] Thus, a facile mechanism through which morphology can be modulated can be realized through parameter selection. Increasing jet impingement velocity leads to a finer distribution of well-defined channels, while material removal increases with greater current densities and lower feed rates evidenced by deeper channels. The latter is intuitive considering the respective increase in the linear charge density (total charge transferred per length of toolpath, Figure 5d), ranging from 0.5-3.1 C mm −1 at the parameters tested. Faster feed rates (ν f > 0.3 mm s −1 ) lead to discontinuous pitting across the toolpath length at linear charge densities < 0.5 C mm −1 , while deterministic grooves, with broadly parallel side walls were observed > 3.1 C mm −1 .
Pitted area per toolpath length increases with linear charge density at each given feed rate, correlating with material removal (Figure 5e; Figure S4, Supporting Information). However, between different feed rates this relationship does not hold due to breakdown propensity increasing with applied current density (Figure 5f). Lower applied currents and potentials stimulate fewer discrete breakdown points and charge is transferred through a smaller number of exposed sites. This promotes the stability of pre-existing pits and suppresses pitted area regardless of linear charge density. areas are suppressed with increasing jet velocity for a given current density and feed rate. Increasing jet velocity decreases the thickness of the diffusion layer and increases the supply of iodide to the anode surface, which suppresses pitted area. The rate of iodide oxidation (Equation 2) and the corresponding limiting current increase, reducing the quantity of material removal. [55] Translating the jet leads to more complex features such as anastomoses where microchannels propagate into one another  . e) For a given feed rate, pitted area increases with linear charge density, while increasing jet velocity suppresses pitted areas. Density of unique features, such as branch points f) and termination points g) decreases with linear charge density as the marks become more deterministic. h) A taxonomy of vascular type features that could be applied to further image recognition. i) Photograph of an engraved physical tag with C 6 toolpath symmetry in the form of a snowflake. j) Optical micrograph (20x objective) of red boxed area in (i). k) Optical micrograph (100x objective) of boxed area in (j). l) Photograph of the snowflake pattern along with a process flow showing an individual part identification routine with an internet-connected smartphone to compare with a library or database of watermarks. Adoption of a blockchain would also allow the creation of a digitalised twin of a given watermark to provide a public record of ownership and authenticity m). Scale bars are 0.5 mm unless stated.
to form a junction. Anastomosis means that there will "not necessarily" be at least n + 1 terminations for each branch point. This is in contrast to patterns marked without translation (Figure 4c), where pits propagate radially outward directed by the jet flow, such that pit divergence is intrinsic to propagation. By translating the electrochemical jet, an additional growth bias can be applied, and the propagation vector can be manipulated. Anastomoses are suppressed as the jet velocity is increased, i.e., the radial component of the propagation vector is increased, for any given current density and nozzle translation velocity. To accommodate the morphological complexity, a feature taxonomy can be proposed; again, we are inspired by the classification of fingerprint minutiae (Figure 5h). These can be exploited to create more complex verification routines.
The apparent disorder is dependent on the length scale of observation. In this sense, the vascular arrays can be considered to have a fractal nature. As an example, micrographs of a 'snowflake' electrochemical pattern marked with a rotationally symmetric toolpath at different magnifications are shown in Figure 5i-k. There is apparent self-similarity as magnification is increased-disorder crosses length scales. This is unsurprising as the flow mechanisms that dictate the pit propagation vector also cross these length scales. The complexity observed at different length scales provides more evidence that the patterns created through this method cannot realistically be duplicated, even via equivalent fabrication routines. Engraved unclonable watermarks like the coniferous branch (Figure 1i), or the 'snowflake' pattern (Figure 5l), can be associated with a single, tangible part or item. This can be matched against a digital equivalent (i.e., a 'twin') within a library through an image recognition routine. This digital equivalent can also be cryptographically encoded, for example as a non-fungible token (NFT) incorporated into an immutable ledger (for instance a public blockchain) that could allow confirmation of authenticity and a record of ownership, or maintenance history, of a physical part within the digital space (Figure 5m).
To be deployable, watermarks must be resistant to duplication from counterfeiters. Watermarks visible to the naked eye are imageable (we consider this to be an asset) and theoretically reproducible by other engraving methods. While this is not a steganographic approach, there are few methods capable of engraving surfaces to this filigree. Shear-based micromachining used in conventional engraving cannot approach these length scales. Thermal-based methods like lasers and high energy electron beams might achieve structures approaching these length scales, but often result in surface damage like burrs and a heat affected surface that would be readily apparent on inspection by the end user. Consider the etched surfaces in Figure 2h-j; these are characteristic of anodic dissolution that is an athermal process in which the underlying microstructure is revealed but not altered. A focused ion beam can achieve the resolution, but over an unrealistic timescale (up to 0.5 µm 3 nC −1 ). [56] In any case, the characteristic etch surface would not be apparent.

Thermal Stability of Patterns
To appraise the thermal stability of the watermarks at elevated temperature in an oxygenated atmosphere, patterns were exposed to heating cycles (100-750 °C) under exposure to standard atmospheric composition (Figure 6a). Figure 6b shows a six-fold symmetric pattern prior to sequential heating cycles, while the micrograph in Figure 6c shows an individual spur of this pattern in the initial condition. Figure 6d-g show the same spur after sequential heating, such that the final micrograph shows this spur at the culmination of all of the preceding heating cycles (Figure 6g). Uniform thermal oxidation (yellowing) of the surface is shown at the 500 °C condition, where the coloring is indicative of a thin oxide film (<1 µm). Surface oxidation is extensive and discontinuous after the 750 °C condition. This oxidation behaviour is typical of this Ti-alloy, which is commonly used in low to medium temperature applications (up to 400 °C) in oxygenated atmospheres.
Features appear dilated at the highest temperature condition tested (750 °C), which may affect pattern readability. This is enlarged in Figure 6h,i for the initial and 750 °C conditions, respectively, where the filigree of the branches is reduced in the latter case. This apparent dilation results from oxidation at the edges of the microchannels, affecting the reflectivity and apparent contrast. This is evidenced by the accompanying secondary electron micrographs of two different regions of the pattern where microscale structuring appears to be lost (Figure 6j,k). In both cases, the morphology of these features is consistent with the virgin condition shown in Figure 6h, despite the severe degradation of the surrounding surface of the material. In any case, the patterns remain visible after even the most extreme thermal cycle, which lies well beyond the operational range of this alloy. This indicates that a pattern may still be readable upon exposure to conditions that would nominally lead to the failure of a given part.
The readability of 15 static patterns (without jet translation) after an aggressive thermal cycle (550 °C, 1 hr) was assessed by imaging before and after heating, where images are shown in both conditions in Figure 6l for a typical pattern. Patterns were marked at similar parameters to those shown in Figure 4 (J = 100 A cm −2 , 0.75 s, ν f = 24.1 m s −1 ). Figure 6m shows these images after registration and binarization, while the correlation coefficient for the polar data strings is shown for all 15 patterns in Figure 6n. The results shown in Figure 6m,n indicate that patterns are identifiable even after exposure to conditions that severely degrade the material. This is intuitive as the etching process is not altering the composition or the microstructure of the material, thus the thermal behaviour (including decomposition and oxidation behaviour) is retained. This is one intrinsic advantage of 'top down' processing methods like etching as compared with methods that coat the part with a discrete material. Likewise, failure mechanisms like thermally induced delamination, driven by differential thermal expansion between the substrate and a pattern fabricated from a different material, are eliminated.

Conclusion
To conclude, we have shown that turbulent and diffusive mass transfer regimes can be exploited in combination with stochastically distributed defects in the surface of robust engineering alloys to allow the electrochemical engraving of dendritic structures under the application of an electrical field. This is achieved with a factory floor compatible machine tool without modification. The watermarks are of the length scale that can readily be imaged with smartphone cameras. While the flowdirected structures are chaotic, they abide by a set of design rules and can be propagated along a deterministic toolpath vector, while retaining individuality. As pattern formation relies on the interplay between the chemical environment adjacent to the anode and the anode material, including the surface oxide, discrete parameters will be material dependent. However, these adaptions are trivial to implement in a machine tool, wherever there is a desire to authenticate a given part. We have shown that logical data strings can be extracted from these patterns using computationally inexpensive feature extraction methods, and these strings can be rapidly cross correlated to discriminate one pattern from the next. As such, we have demonstrated the principle of rapidly direct-writing unclonable marks that can be applied toward authentication of high-value parts and consumer goods that are vulnerable to counterfeiting and we have shown that these marks are identifiable after exposure to thermal conditions that lead to the severe degradation of the material.

Experimental Section
Materials: Titanium alloy sheet (Ti-6Al-4 V, 3 mm) was used as the substrate after sequential grinding and polishing (up to 1 µm diamond) and rinsing in deionized water. Marking was performed using a NaI electrolyte (2.0 m, Fisher Scientific), which was made up using the as-received salt dissolved in deionized water. Electrolytic conductivity was measured to be 154.9 ± 0.1 mS cm −1 at 23.1 ± 0.1 °C, while density was measured volumetrically to be 1.22 g L −1 . Marks were also fabricated using a NaCl electrolyte (Fisher Scientific), of similar electrolytic conductivity (147.7 ± 0.1 mS cm −1 at 21.8 °C).
Fabrication of Vascular Arrays: Marking was performed using a 3-axis machine described in previous work. [57] All experiments were carried out using a 0.51 mm diameter nozzle, at controlled inter-electrode gap (IEG) length (0.50 mm). Striations were processed along 5.0 mm path lengths. Experimental variables are as follows: Current densities, J, I A N       where I is applied current and A N is cross-sectional nozzle area of 50 -150 A cm −2 were generated with applied currents of 0.10-0.31 A. Feed rates, ν f , of 0.1, 0.2, and 0.3 mm s −2 , and jet velocities, ν j , 18.6 ± 0.15, 20.4 ± 0.10, 23.4 ± 0.05 m s −1 were used. In all cases, high Reynolds number (>10000) turbulent flow was expected at the nozzle outlet.
Thermal Cycling: Two thermal cycling experiments were undertaken in this study. First, Sequential heating of one sample was undertaken at over a range of temperatures (100, 250, 500, and 750 °C) for 0.5 hr. The sample was imaged between each cycle. Second, 15 individual patterns were heated at 550 °C for 1 hr. In each case, the samples were placed into a preheated furnace, removed after the allotted exposure time, and allowed to air cool. The atmosphere was not controlled during the thermal cycling treatment and was taken to be the standard atmospheric composition.
Characterization: Surface geometries were characterized using focus variation microscopy (Alicona G4, 50x and 20x objectives) and scanning electron microscopy (Phillips XL-30). Topography datasets were processed using MountainsMap software to extract removal volumes. Smartphone images were acquired with a Huawei P20 Pro retail handset using the in-built camera without modification. Potential traces were acquired at constant sampling rate (10 ms), where the potential was measured directly across the power supply. Potential traces were filtered by a five-point moving average prior to plotting. Final potential (Figure 3d) was measured as the average measured potential of the last 150 ms of machining for interactions >1.00 s).
Circularity data, 4 Area (Perimeter) 2 π       , were extracted from photographs of marks. Only features > 150 µm 2 were included in measurements (Table S1, Supporting Information). These data were weighted for pit surface area, such that the contribution of larger pits was greater for a given mark. This was performed in ImageJ after binarization. Branch and termination points were extracted from binarized photographs using an inbuilt skeletonization morphological function in MATLAB.
Logical data strings were constructed from the radial distances and polar angles of the extracted termination points, where a positive bit indicates one or more termination points at a given sampling interval. The entropy, H, of these strings was calculated according to Equation (5): where P x is the probability of a positive bit on string x. Hamming distances were calculated as part of a logical XOR operation between two data strings. For polar strings that were rotationally dependent, a circular shifting function was applied such that each string was iteratively shifted by one bit against the next string and the minimum Hamming distance was quoted as the result. Pearson's linear correlation coefficient between two strings a and b was calculated according to Equation (6): where x is the mean of a or b. To aid visual perception, [58] intensity figures were colored according to the "viridis" perceptually uniform color map. [59] Statistical Analysis: Histograms in Figure 4b were generated from the radial distances of 3993 termination points and 527 branch points. The area histogram in Figure 4d was generated from the results of 40 individual patterns. Hamming distance histograms in Figure 4j were generated from the polar and radial data strings extracted from 60 individual pits. The associated curves in each histogram were the Weibull distribution of these data. Histograms were generated in Origin Pro 2021. All other data processing was performed in MATLAB 2021.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.