domain: spi spi

text:
Jeroen Wackers
email: jeroen.wackers@ec-nantes.fr

Ganbo Deng

Emmanuel Guilmineau

Alban Leroyer

Patrick Queutey

Jules Richeux

Michel Visonneau

TOWARDS INDUSTRIAL USE OF ANISOTROPIC HEXAHEDRAL MESH ADAPTATION FOR HYDRODYNAMIC FLOWS

Keywords: mesh adaptation, anisotropy, hexahedral, mesh deformation, automatic setup, ship hydrodynamics

This paper discusses the preparation of the flow solver ISIS-CFD for industrial simulation of free-surface hydrodynamic flows with adaptive mesh refinement. To accurately capture the free surface and boundary layers, anisotropic hexahedral mesh refinement and derefinement is performed by local cell division. The adaptation is based on metric criteria which combine Hessians of the flow variables and free-surface detection. Analytical mesh deformation is used to account for body movement, while a dedicated local deformation algorithm handles the projection of the adapted mesh onto the geometry. It is shown that computational setup protocols can be defined to choose the simulation parameters, without user intervention, for large classes of simulations. Tests on different ships and lifting hydrofoils show the efficiency and versatility of the method.

INTRODUCTION

Industrial application of mesh adaptation is gaining momentum, thanks to the increasing power and reliability of mesh adaptation algorithms on one side and the growing industry awareness of their benefits on the other side. Mesh adaptation was identified as a future key technology in NASA's CFD Vision 2030 [START_REF] Slotnick | CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences[END_REF] and is currently used in aerospace by for example NASA and Boeing [START_REF] Balan | Verification of Anisotropic Mesh Adaptation for Turbu-lent Simulations over ONERA M6 Wing[END_REF][START_REF] Park | Verification of Viscous Goal-Based Anisotropic Mesh Adaptation[END_REF], and Safran [START_REF] Mercier | Dynamic adaptation of tetrahedral-based meshes for the simulation of turbulent premixed flames[END_REF][START_REF] Nardoni | Bodyfitted discretization, accurate interface tracking and geometrical constraints in topology optimization[END_REF]. For hydrodynamic simulations, mesh adaptation is used with flow solvers such as OpenFOAM [START_REF] Eskilsson | A mesh adaptive compressible Euler model for the simulation of cavitating flow[END_REF][START_REF] Wang | Numerical investigation of unsteady cloud cavitating flow around the Clark-Y hydrofoil with adaptive mesh refinement using OpenFOAM[END_REF], ReFresco [START_REF] Windt | Adaptive mesh refinement in MARIN's viscous flow solver ReFRESCO: implementation and application to steady flow[END_REF], StarCCM+ [START_REF] Yilmaz | An improved Mesh Adaption and Refinement approach to Cavitation Simulation (MARCS) of propellers[END_REF], and FINE/Marine [START_REF] Hildebrandt | Business and technical adaptivity in marine CFD simulations: Bridging the gap[END_REF][START_REF] Yvin | Tip vortex cavitation inception without a cavitation model[END_REF][START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF].

The ISIS-CFD Navier-Stokes solver, developed at Centrale Nantes / CNRS and commercialised by Cadence Design Systems as part of FINE/Marine, includes an adaptive grid capability since 2009. Industrial use of this feature rose sharply around 2018 when the software had matured and accumulated experience al-lowed the user interface to be simplified; it is now used routinely by about half of the industrial users. The principal application is mesh refinement to capture water surfaces [START_REF] Wackers | Industry application of adaptive grid refinement -the case of FINE/Marine[END_REF].

Despite research progress [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF][START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF][START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF], Hessian-based mesh adaptation with ISIS-CFD has so far seen little industrial application. Routine simulation of realistic flows, where the entire mesh is created using adaptive refinement, presents several mesh-related challenges:

• Most hydrodynamic flow solvers use hexahedral meshes with highly stretched cells to capture boundary layers, wakes, and the free surface. Thus, for mesh adaptation, anisotropic refinement of hexahedral meshes is essential.

• Even in steady simulations, the position of a ship is rarely fixed; at the very least, the freedom to move to an equilibrium position is often needed.

Furthermore, refined meshes must be adjusted to conform exactly to the CAD topology. Both issues require mesh deformation procedures.

• For industrial simulation to be successful, the simulation setup should be simple and user guidelines should be reliable and able to handle a large variety of different cases, without requiring expert knowledge.

This article discusses the work done to prepare Hessian-based adaptation for industrial application.

The mesh refinement and deformation capabilities are available today, while on-going work intends to show that straightforward computational setup protocols can be established for all the typical classes of hydrodynamic simulations. As examples, protocols have been presented for displacement hulls like cargo ships and military vessels [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF] and for lifting hydrofoils used on modern high-speed sailing craft [START_REF] Richeux | Towards the automation of adaptive RANS simulations for hydrofoils[END_REF].

The paper is organised as follows. Section 2 presents the flow solver and its specific features for hydrodynamic problems. Section 3 then describes the grid adaptation method, concentrating on those aspects that are particular for hexahedral meshes and hydrodynamic flows. Mesh deformation algorithms are discussed in section 4, while section 5 studies automatic computation setup and parameter choices for large classes of computations. Finally, in section 6 the method is demonstrated for ships and hydrofoils. The conclusion (section 7) discusses the readiness and the usefulness of the procedure for industrial application.

THE ISIS-CFD FLOW SOLVER

ISIS-CFD is an incompressible unstructured finitevolume solver for multifluid flow [START_REF] Queutey | An Interface Capturing Method for Free-Surface Hydrodynamic Flows[END_REF]. The velocity field is obtained from the momentum conservation equations and the pressure field is extracted from the mass conservation constraint transformed into a pressure equation. These equations are similar to the Rhie and Chow SIMPLE method [START_REF] Rhie | A numerical Study of the Turbulent Flow Past an Isolated Airfoil With Trailing Edge Separation[END_REF], but have been adapted for flows with discontinuous density fields. The method features turbulence models like Menter's k -ω SST [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] but also an anisotropic EASM model [START_REF] Duvigneau | On the role played by turbulence closures in hull shape optimization at model and full scale[END_REF] and hybrid RANS/LES models of the DES family. Imposed and resolved motions of rigid or flexible bodies can be computed as will be shown in section 4, while large-amplitude multibody motion is possible with sliding or overset meshes. The code is fully parallel using the message passing interface (MPI) protocol.

Free-surface flow is simulated with a water-air mixture model [START_REF] Queutey | An Interface Capturing Method for Free-Surface Hydrodynamic Flows[END_REF]. Instead of explicitly reconstructing the water surface, its position is obtained implicitly from a conservation equation for the volume fraction of water. This equation is discretised with specific compressive discretisation schemes to keep the transition between air and water close to a discontinuity (see figure 1). At the surface, regular and fine meshes are required to accurately resolve the volume fraction. Since free-surface flow simulations need regular meshes for good accuracy, ISIS-CFD simulations are mostly performed on unstructured hexahedral meshes like the ones generated by Hexpress from Cadence Design Systems. These meshes feature local refinement with hanging nodes and body-fitted boundary layer grids.

An example of such a mesh is shown in figure 3.

MESH ADAPTATION

Anisotropic hexahedral refinement

The natural mesh adaptation approach for hexahedral meshes is adaptive refinement by local cell division, starting from a coarse initial mesh (figure 3). Hydrodynamic flows contain many flat or stretched flow features, such as vortical wakes or the water surface.

To efficiently capture these, anisotropic mesh adaptation is a necessity [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF]. Therefore, the adaptation in ISIS-CFD [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF][START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF][START_REF] Wackers | Adaptive grid refinement for hydrodynamic flows[END_REF] is performed by anisotropic division: cells can either be split in one or in several directions. For unsteady flows or to accommodate the convergence of a steady flow, existing refinement can be undone. This derefinement is also anisotropic and refinement can be created and removed in a different order [START_REF] Wackers | Adaptivity for complex flows[END_REF].

During a simulation, the refinement procedure is called repeatedly, until the flow has converged and the mesh is no longer changed by the adaptation steps. Since a cell can only be refined or derefined once during each step, the adaptation is an iterative process which requires many calls (40-50) to converge. Therefore, the flow is not fully solved on each adapted mesh. Instead, the adaptation procedure is called every 10 -25 time steps and the flow converges together with the mesh.

Division-based mesh adaptation has its limits. Specifically, the initial grid has a large influence on the adapted mesh, since (1) the division procedure can only reduce the initial cell sizes by a power of two, so not all cell sizes can be obtained, and (2) the adaptation cannot modify the orientation of the cells. This means that the initial grid with its body-aligned boundary layer grid has to be created with care; coarse but orthogonal grids and thick boundary layer meshes produce optimal results. An illustration of an initial and adapted ship mesh is provided in figures 3 and 2, while figure 4 shows a mesh for a hydrofoil.

Refinement criterion

The decision where to refine is based on a metric refinement criterion [START_REF] George | Creation of internal points in Voronoi's type method. Control adaptation[END_REF][START_REF] George | Delaunay Triangulation and Meshing -Application to Finite Elements[END_REF][START_REF] Alauzet | Feature-based and goaloriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace[END_REF], which is applied to divisionbased refinement of hexahedrals as follows. Starting from a symmetric tensor field C(x, y, z) computed from the flow, the grid is adapted until the three dimensions di,j of each hexahedral cell i, which are the vectors between the face centres for the the three pairs (j = 1, 2, 3) of opposing faces, satisfy:

|Ci di,j| = Tr, ∀i, ∀j ∈ (1, 2, 3), (1)
as closely as possible. Here |•| indicates a vector length.

The constant threshold Tr functions as a global specification of the mesh fineness. Cells are refined in a direction j whenever the modified distance |Ci di,j| exceeds Tr and derefined whenever it is below Tr/2.5, to make sure the derefined cells are not immediately refined again. This definition leads to anisotropic refinement whenever the cells are aligned with a non-uniform metric, which is typically the case at the water surface.

Although the formulation (1) differs slightly from the approach of George et al. [START_REF] George | Creation of internal points in Voronoi's type method. Control adaptation[END_REF], who use a metric MP to directly compute a modified vector length u T MP u, this difference is cosmetic. If C is chosen such that C C = MP , the resulting modified distances are the same between the approaches.

To create fine meshes entirely using adaptive refinement, the refinement criterion C must react to all the flow features which are relevant for ship resistance.

For the orbital flow fields in the waves, the Hessian of the pressure is a suitable indicator [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF]. The accurate resolution of the boundary layer and wake is needed to correctly predict both viscous and pressure forces, which suggests refinement based on the velocity Hessians as well. Thus, the criterion is based on the Hessians of both the pressure and the velocity, weighted in the way in which they appear in the flux [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF]. A common weight ρV is assigned to all the velocity Hessians, where

V = √ u 2 + v 2 + w 2 : CH = max H(p) , ρV H(u) , ρV H(v) , ρV H(w) ϕ . (2)
H is the Hessian operator.

• indicates a matrix having the same eigenvectors as the original one and the absolute values of its eigenvalues; the power ϕ is again applied to the eigenvalues, where ϕ = 0.5 is used in practice. Finally, the maximum of two tensors is computed following the procedure in [START_REF] Wackers | Adaptive grid refinement for hydrodynamic flows[END_REF], which is based on [START_REF] George | Delaunay Triangulation and Meshing -Application to Finite Elements[END_REF]. The convection-derived weight ρV for the velocity Hessians is open to discussion, since these mainly indicate diffusion-based errors. This is a subject of on-going study.

The power ϕ can be used to give more or less weight to the strongest flow features. For ϕ = 1 2 , the equivalent metric of [START_REF] George | Creation of internal points in Voronoi's type method. Control adaptation[END_REF] is [START_REF] George | Creation of internal points in Voronoi's type method. Control adaptation[END_REF] themselves use. Hessian metrics are linked with interpolation error control and, following [START_REF] Alauzet | Feature-based and goaloriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace[END_REF], the use of a Hessian without weighting by its determinant implies a minimisation of the interpolation error in the L ∞ norm, i.e. an equidistribution of the error. However, for our hexahedral finite-volume discretisation, the interpolation error may not be proportional to the discretisation error, so the latter is not equidistributed.

MP = C 2 H ∼ H(p) , H(u) , which is what
Finally, a smooth anisotropic mesh at the free surface is required for the volume fraction equation. Since the pressure gradient is discontinuous at the free surface due to hydrostatic forces, the pressure Hessian is undefined there. To solve this issue, the Hessian criterion is extrapolated from below through the surface region and a second criterion [START_REF] Wackers | Adaptive grid refinement for hydrodynamic flows[END_REF] creates the refined mesh at the free surface: this criterion is non-zero whenever the water volume fraction α is neither 0 nor 1. Directional refinement normal to the surface is obtained from normal vectors v = ∇αA/|∇αA|, where αA is a smeared volume fraction field created with Laplacian smoothing. The criterion then becomes:

CS = v ⊗ v if 0.1 ≤ αA ≤ 0.9, 0 otherwise. ( 3
)
The choice is based on αA to create a buffer layer of a few cells around the surface position for safety.

The two tensor criteria are combined into one [START_REF] Wackers | Combined refinement criteria for anisotropic grid refinement in free-surface flow simulation[END_REF] by taking a weighted maximum of the tensors. Since the free-surface criterion always has a unit eigenvalue, a weighting constant c is applied only to the Hessian criterion:

CC = max (CS, c CH ) . (4)
In practice, choosing c directly is inconvenient, so the user specifies separate thresholds TrS and TrH for the free-surface and Hessian criterion. These are then used to set Tr = TrS and c = T rS T rH . The choice of TrS and TrH is discussed in section 5.2.

In figure 2, combined refinement at and below the surface can be seen in waves, while the mesh around the hull is refined at the bow and in the boundary layer.

Figure 4 shows refinement in the pressure peaks at the leading and trailing edge, as well a anisotropic refinement in the wake.

Quality measures

Certain measures are applied in the mesh adaptation process, to preserve the quality of the adapted mesh [START_REF] Wackers | Adaptive grid refinement for hydrodynamic flows[END_REF]. The two main aspects for this article concern post-treatment of the refinement criterion. First, the eigenvalues of the criteria C can be limited to a userspecified value. This imposes a minimum cell size, below which the cells are never refined. This option is useful, for example, to prevent spurious refinement near flow discontinuities. Furthermore, the metric can be set to zero in the x-, y-and z-direction on the outside of limiting boxes, to prevent refinement close to the border of the domain. The specific use of these features is explained in section 5.

Other quality measures act directly on the refined mesh, like the difference of one refinement level allowed between neighbour cells, or the obligation to keep the same refinement in a column of boundary layer cells. For brevity, these are not further described here.

MESH DEFORMATION

For realistic flow simulation in hydrodynamics, mesh deformation is essential. The first reason for this is, that the forces on a ship are highly dependent on its position with respect to the water surface: a ship that sits slightly deeper in the water has a much higher drag. Thus, even for steady simulations, it is more realistic to let a ship move to its dynamic equilibrium position, than to fix the position. For fast ships and hydrofoils, this can lead to significant motion. Furthermore, for hydrofoils, naval architects are observing that the deformation of the geometry needs to be taken into account to get accurate forces. This requires mesh deformation to follow the motion of rigid and flexible bodies. Finally, division-based mesh refinement creates meshes where the newly inserted nodes are not body-fitted. To correct this, the mesh is locally deformed to project it back onto the geometry.

Since these mesh deformations do not have the same requirements, different deformation methods are employed in ISIS-CFD for rigid-body motion, flexiblebody deformation, and for the projection of adapted meshes. This section briefly describes these methods. For the coupling with the flow equations, ISIS-CFD uses an Arbitrary Lagrangian Eulerian formulation [START_REF] Queutey | An Interface Capturing Method for Free-Surface Hydrodynamic Flows[END_REF] which allows to solve the flow on moving meshes.

Rigid-body motion

For rigid-body ship motions [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fishlike body[END_REF], the challenge is to preserve a good mesh at the free surface. Therefore, rigid-body motions may be divided in those that do not displace the free surface (x-and y-translation, z-rotation, with z the vertical axis), and those that do (z-translation, x-and y-rotation). For single-body flows, the former can be treated by rigid-body displacement of the entire mesh. The latter require mesh deformation where the outer boundaries of the mesh are not deformed, to make sure the mesh at the free surface moves as little as possible. If the free-surface grid is created with adaptive refinement, z-translation can also be treated with rigid-body movement, but xand y-rotation always require deformation.

The deformation is performed with analytical weighted regridding: the required motion is applied to all the nodes of the mesh, multiplied by a weighting coefficient w(x) which goes to zero on the outer boundaries so that only the inner part of the mesh moves (figure 5). An initial coefficient w0(x) is computed as a solution of Laplace's equation with Dirichlet boundary conditions 1 on the body surface and 0 on the outer boundaries. This coefficient is then modified to produce the actual weighting: w(x) = max (w0(x) p /cmax, 1) .

(

) 5
The power p reduces the zone where the mesh deforms significantly, to leave the free surface far away from the ship undisturbed. The factor cmax on the other hand, imposes a zone of rigid-body motion w = 1 close to the hull, to preserve the shape of the boundary layer mesh. For ship simulations with mesh refinement, the default values are p = 3 and cmax = 0.8.

Flexible-body deformation

The deformation of flexible bodies cannot be performed with the analytical weighting presented above.

To perform fluid-structure interaction (FSI), two deformation approaches are included in ISIS-CFD. The first, presented by [START_REF] Durand | Light and flexible Fluid/Structure Interaction, application to sailing boats[END_REF][START_REF] Durand | FSI investigation on stability of downwind sails with an automatic dynamic trimming[END_REF], is to divide the mesh into layers around the hull. The displacement of each face in a layer is copied to the layer around it, defined using the cell-face connectivity, with a diffusion effect to homogenise the deformation and a weighting similar to the one above, to obtain zero deformation on the outer boundaries. The second approach, used with the modal FSI mentioned in section 6.4, is to deform the mesh with radial basis functions distributed over the geometry, to distribute the motion of each of these parts [START_REF] Mouton | Towards unsteady approach for future flutter calculations[END_REF]. These techniques go beyond the scope of this article and are not described in detail.

Post-refinement projection

In division-based refinement, new nodes are inserted in the centres of existing edges and faces. For curved geometries, these nodes are not placed on the geometry, so the refined mesh has to be deformed to fit the geometry. These deformations are highly local, since each misplaced boundary node is surrounded by existing nodes which are placed correctly on the geometry. The layer-based flexible-body mesh deformation mentioned above can be (and has been) used for this projection, but it is not ideal. First, it is inefficient to use such an expensive global deformation to solve essentially local problems. And second, the method does not propagate the local deformations very far from the boundary nodes, which leads to strongly deformed and potentially twisted cells.

Therefore, a dedicated local mesh deformation has been introduced for post-refinement projection. This method starts with the projection of the boundary nodes onto the CAD geometry, represented by a triangulation. Each node is projected orthogonally onto the triangle which gives the smallest displacement from the current node position. In the second step, the interior mesh is deformed. For this, ellipsoids (figure 6) are defined around each boundary node to project, with their long axes oriented along the normal direction to the geometry and with two short axes, of equal length, parallel to the wall. The ellipsoids are dimensioned such that they contain no other boundary nodes. This is obtained by starting from an infinite ellipsoid for each point and then searching the cells nearby. Each time a boundary node is found, either the long or the short axis length is reduced to fit the ellipsoid surface through this node. The axis to adjust is chosen as the one which maintains the largest volume for the ellipsoid.

Once all the ellipsoids are defined, the boundary node deformations (dx1 in figure 6) are applied to all the nodes in each ellipsoid, with a linear weighting so that the outer boundary of the ellipsoid does not deform. Furthermore, a secondary shear translation dx2 is applied to the interior nodes with a quadratic weighting, such that the interior node closest to the boundary node is projected towards the normal vector. Without this secondary lateral deformation, when faces are not aligned with the wall-normal direction, large normal displacements may reduce the orthogonality of the cells to the point that faces can become inverted. Let r(x) be the non-dimensional radius of a point in the ellipsoid, such that r = 0 in the centre and r = 1 on the boundary. Then the total displacement of x is:

dx(x) = dx1 + dx2r(x)/r2 (1 -r(x)), (6)
with r2 the radius of the closest node.

REFINEMENT PROTOCOLS AND COMPUTATIONAL SETUP

To use these algorithms for routine simulations, default values for the simulation parameters are needed, which must be straightforward, reliable, and valid over a large range of ship lengths, velocities, hull shapes, etc. These guidelines are based on physical arguments, notably dimensional analysis, to reduce the number of input variables that have to be considered, combined with experience to find sensible default settings for the remaining parameters [START_REF] Wackers | Adaptivity for complex flows[END_REF].

Simulation protocols have been defined in [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF] for displacement ships and in [START_REF] Richeux | Towards the automation of adaptive RANS simulations for hydrofoils[END_REF][START_REF] Leroyer | Starting off the right foot with foil sock approach and AGR criterion[END_REF] for lifting hydrofoils. Without repeating these protocols fully, we provide examples here of how the guidelines are defined, focusing on the mesh adaptation. Figures 2 and4 show meshes created with these protocols.

Reference length

For non-dimensional guidelines, a key issue is the choice of the reference length, since many parameters such as the threshold Tr of equation ( 1), but also the flow domain size, the boundary layer mesh, certain motion laws, etc. depend on a length scale. To obtain similar meshes for similar flow features, it is important that the length scale represents the characteristic size of the flow, which is typically the length over which the boundary layer develops. For ships, this is the hull length L which is also the main dimension of the ship. For hydrofoils however, the largest dimension is the span but the characteristic boundary layer length is the chord. Therefore, contrary to existing practices which take the span, [START_REF] Richeux | Towards the automation of adaptive RANS simulations for hydrofoils[END_REF] chooses the maximum chord c as the reference.

Below, all lengths are expressed in terms of L for ships and of c for hydrofoils, to underline that the physics of these cases is different. However, both parameters denote the reference length for their respective cases.

Thresholds

The mesh density is specified here by the threshold Tr.

Many authors adjust such mesh density parameters to obtain a specified number of cells, for example by normalising the metric in such a way that the threshold directly specifies the target number of cells [START_REF] Alauzet | Feature-based and goaloriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace[END_REF].

To capture free-surface waves, we prefer a different approach. For non-adapted ship meshes, established ISIS-CFD guidelines specify target cell sizes at the free surface in terms of the ship length L. To obtain the same behaviour for adaptive refinement, the freesurface refinement criterion is defined with unit vectors in equation ( 3), so its threshold TrS equals the desired cell size normal to the free surface, which is chosen as L/1000 like for non-adapted meshes. For hydrofoils, [START_REF] Richeux | Towards the automation of adaptive RANS simulations for hydrofoils[END_REF] suggests c/64 which gives accurate forces, without reproducing all the details of the free-surface flow.

The cell sizes of the initial mesh (section 3.1) should be chosen such that this target cell size can actually be obtained by cell division.

To keep the criteria compatible, a similar behaviour is sought for the Hessian refinement criterion: the mesh The blue lines show the uncertainty interval of the experimental results [START_REF] Kim | Measurement of flows around modern commercial ship models[END_REF].

density on the geometry should be proportional to the threshold value, independent of the case parameters. Therefore, the criterion is non-dimensionalised using the reference length and velocity. As a result, the choice for TrH no longer depends on the hull length and velocity; only a (weak) dependence on F r, Re, and the hull shape remains, which can be ignored for a specific class of geometries (such as displacement ships, fast ships, or hydrofoils). For ships, the tests in [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF] suggest the range TrH ∈ [0.2L, 0.025L] for coarse to fine grids. For hydrofoils, [START_REF] Richeux | Towards the automation of adaptive RANS simulations for hydrofoils[END_REF] recommends TrH ∈ [0.5c, 0.0625c].

Cutoff filter and limiting box

Adaptive refinement based on the combined refinement criterion captures flow details in the entire domain. However, for the accurate evaluation of resistance, not all these details may be required. Thus, the criterion is modified in two ways.

The resolution can be controlled by specifying a minimum cell size, below which cells are no longer refined.

A large minimum cell size works like a cutoff filter: while the adaptive refinement still captures the main flow features, very fine cells to resolve small details are not created. For the KCS test case of section 6.1, figure 9a shows the dependence of the resistance on the minimum cell size. This parameter (chosen small in earlier tests, around 10 -4 L) has a pronounced influence: for a larger minimum cell size, the number of cells is strongly reduced, for only a modest loss of precision. Also, the solutions for all minimum cell sizes converge to the same value when TrH is reduced. Thus, for resistance simulations, a large minimum cell size, i.e. L/1000 like the free-surface threshold, is preferable.

For hydrofoils, the free-surface size was found to be too large so a minimum size of c/256 is recommended.

Furthermore, resolving the wake and the wave field behind the ship may not be needed for resistance computations. Therefore, limiting boxes were tested for the KCS which forbid horizontal refinement from a given distance behind the stern. Preventing refinement aft of 0.3L behind the stern removes the entire far wake field (figure 8). However, compared to a 2L box which preserves most of the waves, the box at 0.3L does not significantly alter the forces while reducing the number of cells by up to 40% (figure 9b). An even shorter box (0.1L behind the stern) changes the converged resistance, so it is not adopted. The box at 0.3L is the best compromise. Similar tests suggest a box at 1c behind the trailing edge for hydrofoils.

TEST CASES

Without attempting a full test of the simulation approach, this section provides four examples of the variety of test cases which it can handle. Where available, the results are compared with experiments and simulations on non-adapted meshes. More complete tests, including further details of the cases presented here, are provided in [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF] for ships and in [START_REF] Richeux | Towards the automation of adaptive RANS simulations for hydrofoils[END_REF] for hydrofoils.

KRISO Container Ship

A modern cargo ship, the KRISO Container Ship (KCS) [START_REF] Kim | Measurement of flows around modern commercial ship models[END_REF], is simulated in model-scale towed condition at Re = 1.257 • 10 7 and F r = 0.260. The k -ω SST turbulence model with a wall law is used. Figure 11 shows the convergence of the resistance with the Hessian threshold TrH , as well as the numerical uncertainty evaluated with the grid-convergence based method of Eça and Hoekstra [START_REF] Eça | A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies[END_REF] (where TrH is used as the measure for the grid size, instead of an averaged cell size). A power-law like convergence is obtained and the uncertainty estimations, which overlap for the finer grids, make sense. Figure 12 shows that the free surface also converges rapidly and that at TrH = 0.05 it is very close to the finest threshold. Thus, varying

TrH is an effective way to obtain grid convergence and the selected range for TrH is appropriate. The results are compared with simulations on Hexpress meshes created with the C-Wizard automatic setup tool, which provides uniform cell sizes over the hull. Also, a combined Hexpress and systematic grid refinement (SGR) approach is shown, where the three coarsest grids are generated by Hexpress, using manual refinement at the bow and stern. The finest grids are created from these coarser grids by refining all cells once, using the adaptive refinement algorithm. Figure 11 compares the three results. The adapted series converges faster to its final value than the other two. A given numerical accuracy is therefore obtained on coarser grids: the same distance from the converged value is typically obtained with 1.5 to 2 times fewer cells than the C-Wizard series. Furthermore, the adapted and SGR series converge to roughly the same values for the resistance, which underlines the reliability of both series. The meshes of the C-Wizard series are not similar: in Hexpress, the boundary layer meshes become thinner as the outer mesh is refined. Therefore, the series converges to a lower value than the other two.

DTMB 5415

The protocol (developed for the KCS) is applied unchanged to the DTMB 5415 model in towed condition at F r = 0.41 and Re = 1.74 • 10 7 following the measurements of [START_REF] Olivieri | Towing tank experiments of resistance, sinkage and trim, boundary layer, wake, and free surface flow around a naval combatant INSEAN 2340 model[END_REF]. This case represents a slender military vessel at full speed, where heavy wave breaking occurs both at the bow and the stern (figure 13).

The adapted-grid resistance is again compared with C-Wizard uniform grids and meshes with manual refinement and SGR (figure 14). The adapted series converges quickly and produces good uncertainty estimations; however, the converged resistance is smaller than for the experimental results, likely due to the difficulty of modelling the breaking wave behaviour. The non-adapted series seem to agree better with the experiments, but they converge to the same underestimated solution -albeit more slowly. Thus, the faster convergence and smaller uncertainty mean that the adapted-mesh results are the most reliable.

KCS in self-propulsion

To test the versatility of the ship refinement protocol, it is applied to a case for which it was never intended: the KCS in self-propulsion condition [START_REF] Hino | CFD Workshop Tokyo[END_REF] with a bodyforce model to represent the propeller. This is a challenge since the mesh refinement has never been tested for a force field before. Also, a ship propeller is tiny with respect to the ship itself. Can something so small be captured by mesh adaptation that is configured for the flow around the ship? The initial Hexpress mesh is equal to the one from section 6.1; no particular refinement is applied around the actuator disk. Still, the adapted grid (figure 15) captures the propeller flow well, even at the coarse threshold TrH = 0.1L. Thus, refinement based on the ship length scale is appropriate for a detail like the actuator disk. However, the mesh size throughout the disk region is equal to the minimum cell size. While this mesh is fine enough for the actuator disk, it is possible that the minimum cell size prevents the convergence to a mesh-independent solution. All resistance results (figure 16) agree well with the experiments; the difference between the series is similar to the DTMB 5415 for example. These force results confirm that the mesh adaptation protocol can be used for self-propulsion.

Flexible META hydrofoil

Finally, the hydrofoil protocol of [START_REF] Richeux | Towards the automation of adaptive RANS simulations for hydrofoils[END_REF] is tested. Like for the ships, this protocol has been established for rigid geometries. However, carbon-fiber foils undergo significant deformation under load and a consensus is developing among naval architects that this deformation should be taken into account for the evaluation of the forces. Therefore, the hydrofoil protocol is applied here with fluid-structure interaction, accomplished via a modal approach [START_REF] Mouton | Towards unsteady approach for future flutter calculations[END_REF]: a structural solver is used to determine the deformation eigenmodes of the foil, then the flow solver computes the deformation as a linear combination of the first eigenmodes. Thus, no direct coupling of the structure and fluid solvers is needed.

The test case is a hydrofoil for an ocean-going catamaran, with a chord c = 0.90m and span of about 2.2m, simulated at 10m/s. The foil has never been constructed so there are no experiments available. Since no structure was ever designed, the foil is considered to be made of a homogeneous isotropic material with a Young's modulus of 60GPa. This simplification affects only the results, not the simulation procedure.

The fluid-structure interaction leads to significant deformation, as shown in figure 17 where the shape is compared with the undeformed foil shape. Both have been simulated with the same hydrofoil protocol. The resistance (figure 18) shows smooth convergence and similar behaviour for the rigid and flexible foils, with good uncertainty estimations. The uncertainty on the finest (23M cell) grid is about 1.6%. All deformed cases have about 200k more cells than the equivalent rigid case; this is likely because the deformed mesh is not perfectly aligned with the free surface, so it requires more isotropic cells to capture the surface.

Finally, figure 19 shows the difference in the drag between the flexible and rigid case. Although the convergence of this difference is not monotonic, the prediction is acceptable even on the coarsest grids. The more general conclusion is that the simulations allow to consistently determine the difference in performance between very similar geometries, even on coarse grids. This is crucial for automatic shape optimisation and it demonstrates the reliability of the simulations.

CONCLUSIONS AND DISCUSSION

In preparing the industrial application of combined Hessian-based and free-surface mesh adaptation with ISIS-CFD, the first question is if such computations are possible today. For the adaptation algorithm itself, this is probably true since mesh adaptation is already in industrial use. The mesh deformation is a standard industrial technique in ISIS-CFD and has long since proven itself. This means that the generality and robustness of the user guidelines are the crucial issue. For calm-water resistance of displacement hulls, the tests in [START_REF] Wackers | Adaptive grid refinement for ship resistance computations[END_REF] show that the protocol can be applied unmodified to hull shapes ranging from slender to full, and to speeds from slow steaming to military top speeds. Furthermore, the protocol is applied successfully to body-force propeller models, a case for which it was never intended. Larger variations of the test case, such as fast planing hulls, were not tried. However, the hydrofoil test shows that a very similar protocol is adequate for these radically different geometries. Most calm-water cases are probably close enough to what was tested here, that even if the protocols in section 5 have to be adapted for them, it is likely that suitable protocols exist. Thus, routine mesh adaptation for resistance simulations appears realistic today.

Does the proposed mesh adaptation offer advantages with respect to existing meshing approaches? The tests show that adapted meshes provide excellent grid convergence, often obtaining 1% difference with the converged resistance on grids of about 1M cells, which is much faster than on non-adapted meshes. Furthermore, mesh adaptation makes the simulation more reliable, since it automatically ensures that the mesh is suitable for the flow. While the correspondence of the resistance with experiments is often worse than for non-adapted meshes, the differences were traced to problems with the non-adapted meshes, such as insufficient resolution of breaking waves. Thus, this is not a numerical problem but a modelling issue. Finally, the grid convergence in all cases presented here is smooth, which leads to reliable convergence-based uncertainty estimations; this was also observed by [START_REF] Wackers | Can adaptive grid refinement produce grid-independent solutions for incompressible flows?[END_REF] for mono-fluid flows and by [START_REF] Alauzet | Feature-based and goaloriented anisotropic mesh adaptation for RANS applications in aeronautics and aerospace[END_REF] for wings and aircraft using compressible flow and tetrahedral mesh adaptation. Thus, the good uncertainty estimations are likely a property of the metric-based refinement. This is an important advantage for unstructured grids.

Open questions remain and the simulation protocol will probably be improved and extended. Still, mesh adaptation has proved its generality and the potential advantages are obvious, so it appears ready for industrial use. Furthermore, hydrofoil simulation in particular has the advantage of being a relatively new application. Many naval architects do not have established simulation practices and may therefore be more likely to adopt mesh adaptation as a standard practice. Thus, the hydrofoil may be the break-through application which will launch Hessian-based adaptation for hydrodynamic simulation. The future will tell.

Figure 1 :

 1 Figure 1: Water volume fraction on the hull of a ship and its centre symmetry plane. Red: volume fraction 1, blue: volume fraction 0.The unstructured discretisation of all the flow equations is face-based. While the unknown state variables are cell-centered, the fluxes which form the systems of equations are constructed face by face, so cells with an arbitrary number of arbitrarily-shaped constitutive faces are accepted. This enables for example adaptive mesh refinement by local cell division: a cell with a divided face and hanging nodes is treated naturally as several independent faces.

Figure 2 :

 2 Figure 2: An example of an adapted mesh. Top to bottom: free surface, details of the Y -symmetry plane, X-cut at the stern. The case is the DTMB 5415 at F r = 0.41 of section 6.2, with Tr,H = 0.07.

Figure 3 :

 3 Figure 3: Initial mesh corresponding to the bottom image of figure 2.

Figure 4 :

 4 Figure 4: An adapted mesh for a lifting hydrofoil geometry, coloured with the water volume fraction.

Figure 5 :

 5 Figure 5: Deformed mesh for rotation of the central body (left) and the weighting coefficient w(x) (right).

Figure 6 :

 6 Figure 6: Ellipsoid around a boundary node, with wall-normal vector, projection displacement dx1, and orthogonality-preserving displacement dx2.

Figure 7

 7 Figure 7 gives an example of the process. It shows the undeformed mesh with the ellipsoids and the final deformed mesh. Around the leading edge of the airfoil, secondary displacement is visible.

Figure 7 :

 7 Figure 7: Example of a refined grid with ellipsoids (left) and corresponding projected grid (right).

Figure 8 :

 8 Figure 8: KCS, Free-surface elevation with limiting box 0.3L behind the stern (left) and 2L behind the stern (right). Hessian threshold TrH = 0.05L. The dashed line indicates the aft edge of the limiting box.

Figure 9 :

 9 Figure 9: KCS, effect of the minimum cell size (a) and limiting box (b) on the resistance convergence with TrH .The blue lines show the uncertainty interval of the experimental results[START_REF] Kim | Measurement of flows around modern commercial ship models[END_REF].

Figure 10 :

 10 Figure 10: Hull shapes of the KCS (a) and DTMB 5415 (b) models.

Figure 11 :

 11 Figure 11: KCS resistance, grid convergence as a function of the number of cells, with estimated numerical uncertainty. The blue lines show the uncertainty interval of the experimental results [32].

Figure 12 :

 12 Figure 12: KCS, mesh convergence of the free-surface elevation. Top: result for TrH = 0.025. Bottom: the difference of the results for TrH = 0.2L, 0.1L, 0.05L with respect to 0.025L.

Figure 13 :

 13 Figure 13: DTMB 5415, wave pattern for TrH = 0.025L.

Figure 14 :

 14 Figure 14: DTMB 5415, grid convergence and uncertainties for three series of meshes. The blue line represents the experimental result from [34].

Figure 15 :

 15 Figure 15: KCS self-propulsion: axial velocity in the ysymmetry plane for TrH = 0.1L and location of the body force (pink box).

Figure 16 :

 16 Figure 16: KCS self-propulsion: grid convergence of the resistance; the blue line represents the experimental result from [35].

Figure 17 :

 17 Figure 17: Deformed (red) and rigid (green) shape of the META hydrofoil.

Figure 18 :

 18 Figure 18: Convergence of the resistance for the rigid and flexible hydrofoils, with uncertainty.

Figure 19 :

 19 Figure 19: Convergence of the difference in resistance between the rigid and flexible hydrofoil.

ACKNOWLEDGEMENTS

We thank the Institute Carnot MERS (ORUP project) and the research directorate of Centrale Nantes for funding the hydrofoil investigation.


Alice Balfourier

Vladimir Mulens-Arias

Florence Gazeau

Florent Carn

Rational Design of Fractal Gold Nanosphere Assemblies with Optimized Photothermal Conversion Using a Quantitative Structure Property Relationship (QSPR) Approach

Assemblies of plasmonic nanoparticles have been proposed for various applications, including photothermal therapy, exploiting surface plasmon coupling phenomena. However, the rational design of fractal nanoparticleassemblyremainschallenging due to the lack of structural characterizations and modelization of real systems. Here we used the quantitative structure property relationship (QSPR) approach, driven by experimental data and statistical analysis, to establish a relationship between structural descriptorsof fractal gold nanoparticle(GNP) aggregatesand their light-to-heat conversion. A total of 160 assemblies of various size spherical GNPs with different polyelectrolyte chains were synthesized, which differ in their global charge, size, mass fractal dimension, and plasmonic properties. Fifteen independent descriptors of structure and properties were extracted and further analyzed by QSPR. Principal component analysis and multilinear regression reveal that light-to-heat conversion ismainly governed by thestructureof theaggregatesand not by thecharacteristicsof the building blocks. This highlights the key role of the fractal dimension of the aggregate and of the ratio of GNP/polyelectrolyte mass to optimize photothermal effects. Rational criteria to optimize light-to-heat conversion within nonideal fractal assemblies of GNP were identified, relaxingon the choice of other parameters, such as GNP or aggregate size, that can beadapted to the desired biomedical applications.

INTRODUCTION

Nanoparticle assemblies have emerged in the last years as a robust and promising strategy to modulate nanoparticle physical properties. [START_REF] Grzelczak | Directed self-assembly of nanoparticles[END_REF][START_REF] Nie | Properties and emerging applications of self-assembled structures made from inorganic nanoparticles[END_REF][START_REF] Gong | Self-assembly of noble metal nanocrystals: Fabrication, optical property, and application[END_REF][START_REF] Wu | Gold nanoparticles in biological optical imaging[END_REF][START_REF] Pileni | Impact of the metallic crystalline structure on the properties of nanocrystals and their mesoscopic assemblies[END_REF] The emerging collective behaviors of self-assemblies have already proved their potential in the biomedical fields with applications for biosensing, [START_REF] Zeng | A review on functionalized gold nanoparticles for biosensing applications[END_REF][START_REF] Saha | Gold nanoparticles in chemical and biological sensing[END_REF] surface enhanced Raman spectroscopy (SERS), [START_REF] Majumdar | DNAmediated wirelike clusters of silver nanoparticles: an ultrasensitive SERS substrate[END_REF][START_REF] Szekeres | SERS probing of proteins in gold nanoparticle agglomerates[END_REF] drug delivery, [START_REF] Chou | DNA assembly of nanoparticle superstructures for controlled biological delivery and elimination[END_REF] magnetic resonance imaging (MRI), [START_REF] Berret | Controlled clusteringof superparamagnetic nanoparticles using block copolymers: design of new contrast agents for magnetic resonance imaging[END_REF][START_REF] Corr | Linear assemblies of magnetic nanoparticles as MRI contrast agents[END_REF] or photothermal therapy (PTT). [START_REF] Li | Mixed-charge self-assembled monolayers as a facile method to design pH-induced aggregation of largegold nanoparticlesfor near-infrared photothermal cancer therapy[END_REF][START_REF] Kang | Mesenchymal stem cells aggregate and deliver gold nanoparticles to tumors for photothermal therapy[END_REF][START_REF] Sun | Salt-induced aggregation of gold nanoparticles for photoacoustic imagingand photothermal therapy of cancer[END_REF][START_REF] Song | Smart gold nanoparticles for photoacoustic imaging: an imaging contrast agent responsive to the cancer microenvironment and signal amplification via pH-induced aggregation[END_REF][START_REF] Mulens-Arias | Polyethyleneimine-assisted one-pot synthesis of quasi-fractal plasmonic gold nanocomposites as a photothermal theranostic agent[END_REF] Nanoparticle assemblies are most often directed by DLVO (Derjaguin-Landau-Verwey-Overbeek) interactions and embrace a wide variety of mono-or polydisperse, isotropic or anisotropic structures, that could be highly organized or randomly arranged. [START_REF] Bouju | Nonisotropic Self-Assemblyof Nanoparticles: From Compact Packing to Functional Aggregates[END_REF] However, despite the great versatility of nanoparticle assemblies, such systems are still hard to optimize for a desired application. Indeed, the multiparametric relationship between a complex architecture and one resultingphysical property can hardly be predicted by analytical approaches, which makes it difficult to optimize assembly structures. [START_REF] Torquato | Inverse optimization techniques for targeted selfassembly[END_REF][START_REF] Adleman | Combinatorial optimization problems in self-assembly[END_REF] In many fields, the rationalization of structure-property relationships for complex objects has been a huge step to predict new efficient objects, for example, in drug design. [START_REF] Mandal | Rational drug design[END_REF][START_REF] Huggins | Rational approaches to improving selectivity in drug design[END_REF] One of these rational approaches, called quantitative structure property relationship (QSPR), or its biological variant quantitative structure activity relationship (QSAR), is based on statistical modeling and aims to build a mathematical relationship between thefeaturesof diverseobjects(molecules, polymers, crystals, nanoparticles, fluids) and a biological activity (interactionswith asubstrate, cytotoxicity), achemical reactivity, or a physical property. [START_REF] Katritzky | QSPR: the correlation and quantitative prediction of chemical and physical properties from structure[END_REF][START_REF] Schultz | Quantitative structure-activity relationships (QSARs) in toxicology: ahistorical perspective[END_REF][START_REF] Roy | A primer on QSAR/QSPR modeling: fundamental concepts[END_REF] This methodology presents the advantage to be driven by experimental data and computer calculation, limiting the bias due to human analysis or the assumptions inherent to theoretical approaches. On the other side, it requires robust and unbiased descriptors and a largenumber of samples. Such computational approacheshave the potential to improve fit-for-purpose assemblies of nanoparticles, as it can address multiparametric complex questions.

Herewechoseto challengethisstrategyfor arational design of gold nanoparticle (GNP) aggregates optimized for photo-thermal therapy. When closely assembled, thelocalized surface plasmon resonances (LSPRs) of adjacent GNPs are coupled, which induces a longitudinal plasmon at higher wavelengths than the initial transversal plasmon. This longitudinal plasmon is of great importance because it could reach the first tissue transparency window in the near-infrared range (NIR1, 650-900 nm). [START_REF] Ghosh | Interparticlecouplingeffect on thesurface plasmon resonanceof gold nanoparticles:fromtheoryto applications[END_REF] Hence, plasmon coupling could enable usage of GNPs of small sizes in diverse biomedical applications and activating them in the NIR1 window. GNP aggregates have been proposed for photothermal therapy or photoacoustic imaging with a proven efficiency but without any control or study of the aggregate structure. [START_REF] Li | Mixed-charge self-assembled monolayers as a facile method to design pH-induced aggregation of largegold nanoparticlesfor near-infrared photothermal cancer therapy[END_REF][START_REF] Kang | Mesenchymal stem cells aggregate and deliver gold nanoparticles to tumors for photothermal therapy[END_REF][START_REF] Sun | Salt-induced aggregation of gold nanoparticles for photoacoustic imagingand photothermal therapy of cancer[END_REF][START_REF] Song | Smart gold nanoparticles for photoacoustic imaging: an imaging contrast agent responsive to the cancer microenvironment and signal amplification via pH-induced aggregation[END_REF] The relationship between the GNP aggregate structure and their plasmonic properties has been studied theoretically using an analytic approach. Khlebtsov et al. highlighted the importance of the interparticle distance, the presence of short linear chain fragments, and the number of particles embedded in the aggregate. [START_REF] Khlebtsov | Optical amplification of photothermal therapy with gold nanoparticles and nanoclusters[END_REF] However, the consistency between theoretical and experimental resultshasbeen assessed on ideal and monodisperse nanometric aggregates that could not be synthesized in water with biocompatible components and with a scalable protocol. On the contrary, nonideal aggregates have not been subjected to comparison between analytic prediction and experimental observation.

Nanoparticle aggregation also occurred upon cellular internalization of GNPs: GNP sequestration in the cell lysosomes modifies their initial optical properties, increasing or decreasing their therapeutic effects. [START_REF] Espinosa | Cancer cell internalization of gold nanostars impacts their photothermal efficiency in vitro and in vivo:toward aplasmonic thermal fingerprint in tumoral environment[END_REF] Understanding the structure-heating property relationship of complex GNP aggregatesisthusof fundamental interest to design therapeutic agents that do not lose their potential in biological media or after sequestration in cells.

Here we propose a QSPR regression model to establish a quantitativerelationship between the experimentally accessible descriptors of plasmonic aggregate structure and their optical and photothermal properties. Nanoparticle aggregation was trigerred by electrostatic complexation between citrate-coated GNPs (<30 nm) and oppositely charged long polyelectrolyte chains. This versatile approach enables GNP assemblies to be obtained in a wide range of size, mass fractal dimension, and inner structure by a simple mixing strategy. [START_REF] Shi | Shape-Tailored Colloidal Molecules Obtained by Self-Assembly of Model Gold Nanoparticles with Flexible Polyelectrolyte[END_REF] GNP aggregates were characterized by dynamic light scattering (DLS), laser Doppler electrophoresis (LDE), small-angle X-ray scattering (SAXS), cryo-transmission electron microscopy (cryo-TEM), UV-visible spectroscopy, and photothermal measurements, providing a set of descriptors suchasthesize,charge,fractal structure,interparticledistances, and optical and heatingpropertiesof theaggregatesat different times after synthesis. In total, 18 descriptors were extracted from these measurements and gathered in a single database compatible with QSPR tools, including data visualization by principal component analysis and statistical model construction (Scheme 1). These combined experimental and computational approaches allowed us to conclude on an optimized configuration of GNP aggregates for photothermal therapy. It gives an example on how QSPR methodology can be useful to optimize nanoparticle assemblies to optimize a given physicochemical property.

RESULTS AND DISCUSSION

Construction of a Data Set to Describe Gold Nanoparticles Polyelectrolyte Self-Assembly. Three sizes of citrate-coated, monodisperse, and well-dispersed gold nanospheres were synthesized using Turkevich inspired protocols (D = 3, 13, and 30 nm) (Figure S1 in the Supporting Information). [START_REF] Bastuś | Kinetically controlled seeded growth synthesis of citrate-stabilized gold nanoparticles of up to 200 nm:sizefocusingversusOstwald ripening[END_REF][START_REF] Piella | Size-Controlled Synthesisof Sub-10-nanometer Citrate-Stabilized Gold Nanoparticles and Related Optical Properties[END_REF] Three biocompatible polyelectrolytes (chitosan 300 kDa, chitosan 700 kDa, branched poly(ethylenimine) (PEI) 750 kDa) were diluted in an acidic medium to ensure their protonation ([HCl] = 10 -3 ). Then, GNP aggregates were synthesized by an equivolumic mix of the GNP suspension and the polymer solution. Six concentration ratios were tested for each series, in which the concentration of atomic gold was fixed to 16.4 mg/L and the concentration of polyelectrolyte was varied between 3 × 10 -5 and 10 -2 g/L (i.e., polyelectrolyte/gold mass ratio between 0.0018 and 0.6). Thereafter, the aggregates were characterized 5 min and 24 h after the mixingof GNPsand polyelectrolytes. The long time (i.e., 24 h) was the first studied to ensure the stability of the structure of the objects under scrutiny, but we realized that, regardless of their structural stability, aggregates of gold particles sediment due to their density, even for small aggregation number when the particles are large. We then Scheme 1. Scheme of the Experimental and Computational Approaches Combined in the Present Study carried out asecond seriesof measurementsat short times(i.e., 5 min) in order to have sufficiently concentrated suspensions and to take into account the photothermal properties of the largest structural variety of aggregates. We point out that for theaggregateswiththehighest global chargethesedimentation kinetics are not perceptible before 30 min, whereas the aggregates are formed in less than 1 s. The data obtained at long times have nevertheless been kept to contribute to the wide range of behaviors embraced in this study.

First, the charge and size of the formed aggregates were estimated using their electrophoretic mobility (EM) and DLS relaxation times, respectively. DLS correlograms displayed in some samples a double relaxation process, indicating that two populations of objects can coexist. Accordingly, DLS correlograms were thus fitted with a mono-or biexponential fit, to obtain the relaxation times DLS t1 and DLS t2 (with DLS t2 > DLS t1 ), and their relativeamplitude, DLS A1 and DLS A2 (Figure S2A).Indicativehydrodynamic sizecould beextrapolated from the relaxation time measurements but should be taken with caution due to the presence of quickly sedimenting large aggregates (i.e., R >500 nm).

SAXS was performed to assess the fractal structure of the aggregates. At the higher q values, we always observed the sphere form factor of individual GNP building blocks, corresponding to the TEM size. At the lower q values, the curves show a power law shape, observed on at least one q decade, that enable us to extract the fractal dimension, FD, of the assembly (Figure S2B). The fractal dimension in question isthemassfractal dimension describingthespatial distribution of the mass within the aggregates. For an aggregate composed of particlesof massm 0 and sizer 0 , themassof theaggregatem scales with its radius r, measured from any site within the fractal structure, as ( )

m r m ( ) r r FD 0 0 ∝
. In other words, FD characterizestheramification of thecluster and variesbetween FD ∼ 1 for one-dimensional clusters and FD ∼ 2.9 for threedimensional compact clusters. SAXS did not give us access to the size of the aggregates, as no Guinier regime was observed at low q values. The distance between the GNPs inside the aggregates was investigated bycryo-TEM in somechosen aggregatesthat were closeto thechargeneutralityratio.Theinterparticledistanceis known to be a key parameter to predict the optical properties of GNP aggregates,asplasmon couplingisonlyeffectiveunder really short distances. [START_REF] Khlebtsov | Optical amplification of photothermal therapy with gold nanoparticles and nanoclusters[END_REF][START_REF] Louis | Gold nanoparticles for physics, chemistry and biology[END_REF] In addition to the structural characterization of GNP aggregates, their optical properties were investigated. UVvisible extinction spectra display one or two plasmon bands, correspondingto the longitudinal and transversal plasmon. To isolate the plasmon band features, the spectra were treated to removethebaselineof thesignal, that wascalculated thanksto the spectrum of the single GNPs (Figure S2C). Then, the corrected signal was fitted with one or two Gaussians centered on the mean wavelengths UVvis m1 and UVvis m2 (with UVvis m2 >UVvis m1 ), with a width UVvis w1,w2 and amplitude UVvis A1,A2 , thusprovidingthreeto sixnewdescriptors. Astheinfluenceof scattering effects strongly increases with the size of the aggregates, it could lead to an overestimation of the amplitude descriptors. [START_REF] Huang | Gold nanoparticles: optical propertiesand implementationsin cancer diagnosisand photothermal therapy[END_REF] Finally,photothermal measurementswereperformed usinga 680 nmlaser excitation with apower of 6.3 W•cm -2 at a4 cm distance from the sample surface. The temperature of the sample (V =100 μL) was monitored with an infrared thermal camera during the 5 min of irradiation, until the temperature reached a plateau. To quantify the heating capacity of the sample, the specific absorption rate (SAR) was deduced from the initial slope of temperature increase, T t d d , during the first seconds of irradiation, normalized by the ratio between the gold mass concentration (C m = 16.4 mg•L -1 ) and the heat capacity of water (C p = 4185 J•L -1 •K -1 ; Figure S2D). SAR measurements are a reliable evaluation of heating efficiency because the initial slope of the temperature increase weakly depends on thermal exchange with the environment and thus on the particular experimental setup. Beingalso normalized by the mass of materials, it provides avaluable metric to compare the heating efficiencies of different nanostructures.

All of these measurements were then used as descriptors for the aggregates, in addition to the following synthesis parameters: the GNP diameter (DiamNPs), the average molar mass (MassPoly) and the charge density (ChargePoly) of the polyelectrolyte, the ratio of polyelectrolyte over gold mass concentration (RatioPolyGold), and finally the time after mixing (Time). In total, our database includes 54 types of aggregates observed at two times, plus the three sizes of GNP building blocks: it represents 111 samples described by 18 descriptors,assummarized in Table1.Thisscreeningof alarge number of aggregatesenablesusto establish aphysicochemical relation between series obtained with a given GNP size and polyelectrolyte. Physicochemical Behavior of Gold Nanoparticle Aggregates and Their Properties. Analyses of each series of dataobtained 5 min after mixingreveal acommon behavior when varyingthepolyelectrolyte/gold massratio over 3 orders of magnitude, regardless of the GNP size and the polyelectrolyte nature (Figures S3-S10). Figure 1 displays a case study obtained by mixing 30 nm GNPs and 300 kDa chitosan.

Briefly, theformed GNP assemblieschangefromnegativeto positive global charge (-0.03 to +0.01 mm 2 •V -1 •s -1 ) when increasing the polyelectrolyte/gold concentration ratio and cross the neutrality at concentration ratio between 6 × 10 -3 and 60 × 10 -3 , depending on the series (Figure 1A). Neutral and weakly charged aggregates display the highest size and compactness (highest DLS relaxation time DLS t1 and fractal dimension FD; Figure 1B andC). However, the neutrality point does not systematically yield the maximum value for these descriptors. On the opposite, in a large excess of GNPs or of polymers (mass ratio below 1.8 × 10 -3 or above 180 × 10 -3 ),GNP aggregateshavehigher chargesand their structural features are closer to single GNPs.

Close to the neutrality point, the aggregation is very fast (less than 1 min) and shows a quick color change of the solution. Over time, these aggregates tend to coalesce to form larger aggregates that sediment, and the sample is fully destabilized after 24 h. For these samples, a second DLS relaxation time (DLS t2 ) is observed between 20 and 800 ms, confirming the presence of large aggregates. On the other hand, an excess of charge slows down the aggregation process and the color changes only after several hours, concomitant with an increase of the aggregate size (Figure S11).

Cryo-TEM observations performed close to the neutrality ratio validate that all observed GNPs are embedded into an aggregate and depict a branched structure compatible with SAXS measurements (Figure 1D). Inside the aggregates, interparticle distances are generally below 3 nm, but in some cases (less than 5% ), they can reach up to 15 nm. The observation of a low-contrast bridge between the GNPs, probably due to the polyelectrolyte, suggests that GNPs are enclosed in a single aggregate (Figure S12). However, in the large majority of cases, GNPs are close enough to share plasmon resonance oscillations.

This plasmon coupling is reflected in UV-visible spectraby thedecreasein intensityof theplasmon band of thesingleNPs around 525 nmand the appearance of asecond plasmon band between 580and 800nm(Figure1E).Thelatter indicatesthat adjacent GNPs share longitudinal plasmon resonance, in addition to the transversal plasmon resonance observed at 525 nm. This second plasmon band results in an absorption increase at the laser wavelength (i.e., 680 nm) for all samples, including the ones that have a second plasmon band centered on 580 nm. Due to absorption increase, the resulting heating efficiency,that ismeasured asaSAR,increasesfrom20kW•g -1 for singleGNPsto rangefrom50 to 150 kW•g -1 for neutral or weakly charged aggregates (Figure 1F).

Overall, the measurements performed 5 min after the synthesis describe a variety of aggregates from nanometric to micrometric size (roughly 20 to 2000 nm), with a fractal dimension between 1and 2.8,that present acollectiveplasmon peak in the visible to near-infrared range (up to 770 nm). When observed after 24 h, some samples tend to form larger aggregates with time. A general phenomenological behavior could beobserved for most of theseries, revealingamaximum size, compactness, and heatingefficiency of neutral and weakly charged aggregates, and some differences in the aggregation kinetics as afunction of the aggregate charges. However, some interesting variations could also be observed between series. Some descriptors tend to reach a plateau in some series (Figure 1C and F), while others reach a peak (Figure 1B). However, in other series, descriptors have completely different behaviors (e.g., Figure S6B, C, andF). Hence, dataanalysisby this empirical approach does not provide access to fine effects due to the structural variations between the wide range of samples.

Statistical tools are thus required in order to refine our understanding of the structure-heating property relationship. Data mining of this large data set was performed using principal component analysis in order to identify the source of variance between samples. Then, multilinear regression was used to build a predicting model that can link the heating properties of each sample to its conditions of synthesis and structure. As these computational approaches are not commonly used in the nanotechnology field, a brief mathematical description can be found in the Supporting Information. Data Set Exploration by Principal Component Analysis. The experimental data set was first explored by principal component analysis (PCA). As the descriptors RatioPolyGold, DLS t1 , DLS t2 , and SAR have values on several decades, the logarithmic values were calculated for these descriptors. Redundant descriptors were first identified by calculating Pearson's correlation between each pair of descriptors (Figure S13). Four descriptors (DLS A1 , DLS A2 , UVvis w1 , and UVvis w2 ) were suppressed from the data set beforefurther analysis,astheypresent redundancieswith other descriptors (see the Supporting Information).

To discriminate the different types of samples, PCA was performed using the remaining descriptors as input. The first two dimensions, that represent 67.3%of the variance together, clearly reveal three types of behaviors. The latter were extracted by k-means clustering and were noted clusters 1-3 (Figure 2A). Interestingly, the second cluster is oriented along the second plasmon resonance intensity (UVvis A2 ) and the heating efficiency SAR (Figure 2B). This suggests better heating properties than the two other clusters, which is confirmed by Figure 2C, showing the heating efficiency SAR for the different clusters.

Thefeaturesof thissecond cluster arehardlycorrelated with the synthesis descriptors (GNP diameter, type of polymer, polyelectrolyte/gold concentration ratio, or time; Figure S14), indicating that the aggregates included in the second cluster differ mostlybytheir structure.A systematic comparison of the initial data enables us to associate the clustering to three main distinct structures-property behaviors (Figure S15):

• Cluster 1:Theaggregatesizesand thefractal dimensions are close to those of single GNPs. No second plasmon bands are observed (UVvis m2 and UVvis A2 =0). • Cluster 2: The aggregate sizes are higher than those in cluster 1, but no second relaxation time is observed (DLS t2 = 0). UV-visible spectra show two distinct plasmon bands (UVvis A1 and UVvis A2 ≠ 0). • Cluster 3: The aggregates display two relaxation times (DLS t2 ≠ 0). UV-visible spectra exhibit a low plasmon band between 600 and 800 nm but no plasmon band around 525 nm (UVvis m1 and UVvis A1 =0).

These three behaviors show that there is an optimum state of aggregation (cluster 2) that lies between the single GNP behavior (cluster 1) and aggregatesof micrometric size(cluster 3) to optimize the SAR. However, this observation does not provide numerical criteria for this intermediate state, which requires the construction of a model.

Construction of the Heating Efficiency Predicting Model. We then build apredicting model in order to link the structure and synthesis descriptors to the heating efficiency of GNP aggregates.Theinput datawerelimited to the5mintime point, as no measurement was done on the fractal dimension for the24h timepoint.Thetrainingdataset of themodel thus consists of 54 samples, 7 input descriptors (DiamNPs, MassPoly, ChargePoly, EM, RatioPolyGold, DLS t1 , DF), and the output descriptor SAR.

Statistical modeling relies on a chosen mathematical expression in which parameters are optimized in order to fit a training data set with the highest determination coefficient (R 2 >0.7to validatethemodel).Herewestarted fromasimple multilinear model, that has the advantages to be both robust and easy to interpret. We performed a stepwise selection of descriptors in order to limit the number of descriptors that appears in the model. This approach consists of a successive elimination and addition of descriptorsso that thesuppression of a descriptor dramatically decreases the model quality, or reciprocally theaddition of anewdescriptor doesnot improve the model significantly. Moreprecisely, the model is chosen in order to minimize the Bayesian information criterion (BIC), for which the formula is detailed in the Supporting Information.

Our first multilinear regression shows a low determination coefficient (R train 2 =0.62 before stepwise regression, R train+step 2 =0.57 after stepwise regression), which emphasizes the need for amoreflexiblemodel.Thus,weadd quadratic termsfor the parameters that could take more than three values (RatioPolyGold, EM, DLS t1 , DF). As a consequence, the determination coefficient increases up to R train 2 =0.87 before the descriptor selection and R train+step 2 =0.80 after (Figure 3, orange squares). Descriptor selection assumes that linear and quadratic terms are independent variables, to ensure that the descriptors that appear on both forms will not predominate.

Additional testing, such as the normality of the residual distribution, and their distribution through conditions, confirms the mathematical validity of the model (Figures S16and S17).Asthismodel bridgesboththeadvantagesof the multilinear model and a flexibility adapted to the data set, it is selected as the working model.

Once the model has been built, its robustness and its predictability have been evaluated. Indeed, models can overadapt to the training data set, which decreases its predicability. To avoid this effect, two distinct tests have been performed:first an internal validation of the model using cross validation. Briefly, the parameters of the model are estimated on half of the data set and its predictability is assessed on the other half of the data set. We obtain with this method R CV 2 = 0.67, which proves that the model is robust even on alimited number of samples.Second,to beconsidered as robust and predictive, a statistical model must reach a determination coefficient R test 2 higher than 0.6 on the validation data set. [START_REF] Tropsha | The importance of being earnest: validation is the absolute essential for successful application and interpretation of QSPR models[END_REF] For this purpose, a validation data set was constructed, containing aggregates obtained with two new GNP sizes (diameter of 9 and 20 nm, Figure S18) mixed with either the three polyelectrolytes used for the training data set or a new one, polydiallyldimethylammonium chloride (PDADMAC, Mw = 400 kDa). The mixing was made at the same concentration ratios as the training data set, and the characterizations on size and fractal dimension were made 5 min after mixing. This validation dataset contains 48 samples. We emphasize that the validation data set should display descriptor values in the same range as the training data set: same descriptor distribution, GNP diameter between 3 and 30 nm, polyelectrolyte mass between 300 and 750 kDa, polyelectrolyte charge density between -6.2 × 10 -3 and -24.7 × 10 -3 e•mol•g -1 (Figure S19).

The correlation coefficient calculated on the validation data set was R valid 2 =0.64 (Figure 3, blue triangles) and presents an equal residue distribution though conditions (Figure S20). Thus, our model passes all of the validity checks on the mathematical constraints to model predictive values and could be studied more closely fromaphysicochemical point of view.

Polyelectrolyte/Gold Concentration Ratio and Fractal Dimension of the Aggregates Are the Key Parameters for Specific Absorption Rate Optimization. The descriptors that remain in the model after descriptor selection by the stepwise algorithm are the GNP diameter (DiamNPs), the polymer charge (ChargePoly), the polyelectrolyte on gold concentration ratio (RatioPolyGold and RatioPolyGold 2 ), and the fractal dimension (FD and FD 2 ) (Table S1). To evaluate their respective weight in the model, each term was sequentially suppressed and the resulting determination coefficient was calculated on the validation data set (Figure 4). It reveals that the polyelectrolyte/gold concentration ratio and the fractal dimension play major roles in the model, while the polymer charge and the GNP diameter have a minor weight.

Moreover, the two most impacting parameters appear as a second degree polynomial in the model, that present optimum valuesof DF =2 for thefractal dimension and RatioPolyGold = 25×10 -3 for theconcentration ratio between thetwo building blocks. The validity of this prediction is consistent with the observed evolution of SAR with these two parameters (Figure 5A andB). In addition, the optimized concentration ratio calculated by the model is close to the neutrality ratio for all series (Figures 1 andS3-S10). As the concentration ratio could be related to the interparticle distance, this descriptor could account for this lacking structural parameter. [START_REF] Nguyen | Complexation of a polyelectrolyte with oppositely charged spherical macroions: giant inversion of charge[END_REF] The results concerning the fractal dimension appear to be compatible with a previous study on gold nanoparticle assemblies published by Mulens-Arias et al. [START_REF] Mulens-Arias | Polyethyleneimine-assisted one-pot synthesis of quasi-fractal plasmonic gold nanocomposites as a photothermal theranostic agent[END_REF] The authors show these assemblies, that were synthesized by direct nucleation of gold onto a PEI chain, present an increase of their photothermal properties when the fractal dimension evolvesfrom1 to 1.9 accordingto SAXSmeasurements, which strengthen our results. Furthermore, this study evidences the potential of these assemblies for photothermal therapy and photoacoustic imaginginvitroand invivo, which highlightsthe biomedical potential of such objects.

While the conditions of synthesis and the fractal dimension of the aggregates clearly play a role in the model, on the contrary, the building block properties are minor or absent. The polymer mass does not appear in the model, while the polymer charge and the nanoparticle diameter play minor roles. This means that, as long as the concentration ratio and fractal dimension are optimized, the polymer type and GNP size could be chosen independently. To assess this result, we extract thesampleswhich fulfill thefollowingconditions:1.9≤ DF ≤ 2.1 (Figure 5A, red square) and 6 × 10 -3 ≤ RatioPolyGold ≤ 60 × 10 -3 (Figure 5B, red square). Among the 11 samples that fulfill these conditions, and that present significantly higher heating properties than the other samples (p-value=0.0057), thefour polyelectrolytescould befound as building blocks, as well as four of the five GNP sizes (3, 9, 13, and 18 nm). In addition, the first DLS relaxation time varies from130 to 1400 μs between these optimized samples, which covers a hydrodynamic size range from 60 to 680 nm, indicating that the aggregate size plays only asmall role in the heating properties as long as the previous criteria are met.

Thus, the rationalization of GNP aggregate structure with the QSPR approach enables simple and clear criteria to be evidenced that should be fulfilled to obtain optimal heating properties. This identification of criteria relaxes the constraint on the other parameters that govern the system, such as the aggregate size, or the choice of the GNP size and polyelectrolyte used for the aggregation, which enable saferby-design approaches. These parameters can hence be chosen in order to fit other criteriathat are known to depend on their charge and size, such as the biocompatibility of the aggregate, its capacity to avoid or promote elimination by the kidney or the liver, its ability to penetrate deeply inside the tumor, its internalization by tumor cells, or its elimination from the organismafter therapeutic use. [START_REF] Chou | DNA assembly of nanoparticle superstructures for controlled biological delivery and elimination[END_REF][START_REF] Perrault | Mediating tumor targeting efficiency of nanoparticles through design[END_REF][START_REF] Albanese | Effect of gold nanoparticle aggregation on cell uptake and toxicity[END_REF][START_REF] Huo | Superior penetration and retention behavior of 50 nm gold nanoparticles in tumors[END_REF] More particularly, these results raise the attention for the medical use of small nanoparticles, that appear recently to be more biodegradable than large ones. [START_REF] Balfourier | Unexpected intracellular biodegradation and recrystallization of gold nanoparticles[END_REF] Briefly, the author shows that 4 nm GNPs could bealmost entirely degraded by cellsinvitroin 2 months, while 22 nm GNPs are hardly degraded after 6 months. Overall, QSPR provides asafer-by-optimization approach than can really impact the field of nanoassembly design for a large frame of applications.

CONCLUSION

We propose a computational approach based on QSPR regression to model the heating capacities of complex GNP aggregates. For that purpose, we synthesized 108 assemblies using three sizes of GNPs and three types of oppositely charged biocompatible polymers mixed at different ratios. The resulting aggregates cover a range of size from the nanometer to the micrometer scale and that vary from open to compact aggregates. These aggregates were described by 15 independent descriptors that encompass their average size, fractal dimension, charge, and optical and heating properties. Data mining analysis by the empirical approach evidences a general pattern in series,thoughtendenciescannot fullybedepicted by this approach.

Hence, we resort to statistical tools to refine our understanding of this data set. Principal component analysis evidences that samples could be separated into three clusters that are independent of the building block properties, yet highly correlated with heating properties. This indicates that structural propertiesplayamoreimportant rolethan GNP size and polyelectrolytenature, which wasconfirmed by regression. Our analysis emphasizes the weight of the fractal dimension and the GNP/polyelectrolyte concentration ratio in the prediction of aggregate heating properties.

Overall, we validate the use of the QSPR approach on selfassembly structures, and we provide a predictive model to evaluate aggregate heating properties. We evidence that optimizing a few parameters can be sufficient to design optimized functional aggregates. Thislast point enableshaving a choice on other criteria (e.g., on the nanoparticle size) that are important for biomedical application such as the bionano interactions, pharmacodynamics, tissue penetration, intracellular behavior, long-term fate, and safety of these nanodevices.

EXPERIMENTAL SECTION Chemicals. Tetrachloroauric(III) acid (≥99%trace metal basis), trisodium citrate dihydrate (≥99% ), potassium carbonate (≥99% ), and tannic acid were purchased from Sigma-Aldrich. Chitosan 700 kDa was purchased from Acros Organics, and chitosan 300 kDa, PDADMAC 400 kDa, and PEI 750 kDawere purchased fromSigma-Aldrich. Fresh Milli-Q water has been used for all experiments.

GNP Synthesis. The 9, 13, 20, and 30 nm GNPs were synthesized following a Turkevich inspired protocol proposed by ref 30 based on the synthesis of gold seeds of 10 nm followed by asequential addition triggering seed growth. A 97 mgportion of sodiumcitratewasdissolved in 150 mL of Milli-Q water in a two-necked round-bottom flask and heated to boiling under stirring. After 15 min of boiling, 1 mL of preheated 25 mM gold chloride solution was added. After 10 min of reaction,thereaction waseither stopped (for the13nm GNPs) or lowered to 90 °C (for the 30 nm GNPs). The sample was then diluted to enable a quicker growth of the seeds by extracting 55 mL of the GNP solution and the addition of 53 mL of Milli-Q water and of 2 mL of 60 mM sodium citrate solution. A 1 mL portion of 25 mM gold chloride was then added. A second addition of 1 mL of gold solution was realized after 30 min. After another 30 min, the dilution and addition of reagent steps were repeated to obtain the desired size.

The 3 nmGNPs were synthesized in the presence of tannic acid, followingaprotocol proposed by ref 31. A 97 mgportion of sodium citrate and 20.7 mg of potassium carbonate were added to 150 mL of Milli-Q water and 1 mL of 2.5 mM tannic acid solution in a two-necked flask. After heating at 70 °C under stirring, 1 mL of 25 mM gold chloride solution was added. When the solution turned orange, the reaction was stopped.

As shown in Figure S1, these approaches allowed us to obtain stable dispersions of isotropic and monodisperse GNPs with the following radii measured by TEM: D TEM =3.4 ± 0.9, 9.1 ± 1.4, 13.5 ± 3.2, 20.2 ± 2.4, and 30.6 ± 2.6 nm. These GNPs are mainly dispersed as individual objects, as evidenced by UV-visible spectroscopy, DLS, and cryo-TEM. The GNPs were not further treated beforecomplexation and wereused at the lowest concentration obtained after synthesis([Au] =32.8 mg/L). Aggregate Synthesis. Aggregates have been synthesized by an equivolumic mixing of the aqueous GNP solution and theaqueouspolyelectrolytesolution, whosepH waspreviously adjusted to 3 by addition of hydrochloric acid. No salt was added to fix the ionic strength background. Thus, the ionic strength increases with the polycation concentration at fixed GNP concentration. The electrical conductivity was always below 0.5 mS•cm -1 . The polyelectrolyte concentration was modified to obtain final polymer concentrationsrangingfrom3 × 10 -5 to 10 -2 g•L -1 and a gold concentration of [Au] =16.4 mg•L -1 after theequivolumic mix.Theconcentration ratio was evaluated as the ratio of the mass concentration between polyelectrolyte and gold, that is between 0 (GNPs only) and 0.6.

Characterization Methods. Transmission electron microscopy images of single GNPs were obtained with a Tecnai 12 microscope operating at 80 keV (Imagoseine, Jacques Monod Institute, Paris) after deposition of a droplet of NP solution on a hydrophilized grid.

Aggregate characterizations were performed directly after synthesis (less than 5 min after mixing) and after 24 h. Dynamic light scattering (DLS) and laser Doppler electrophoresis (LDE) were performed on a Zeta Sizer NanoZS apparatus (Malvern Instruments), operating with a 632.8 nm laser with P ± 4mW at afixed scatteringangleof 173°at room temperature (25 °C). Size analysis was carried out by converting the measured intensity autocorrelation function into the scattered electric field autocorrelation function using the Siegert relation. Charge analysis was carried out with an electrical field of 25 V•cm -1 oscillating successively at 20 and 0.7 Hz to reduce the electroosmosis effect due to the surface charge of the capillary cell. The particle velocity was measured by LASER Doppler velocimetry.

Small-angleX-ray scatteringexperimentswereperformed on the SWING beamline (SOLEIL synchrotron at St Aubin, France) with the following configuration: sample-to-detector distanceD =6.52 m;beamenergy E =12 keV;λ =1 Å;and q range from 0.0012 to 0.228 Å -1 . The samples were contained in cylindrical glass capillaries of calibrated diameter. The scattered signal was recorded by an Eiger 4 M detector (DectrisLtd, Switzerland) with a pixel size of 75 μm. Preliminary data treatment (angular averaging and normalization) was done using the software Foxtrot developed at the beamline.

GNP aggregates were characterized by cryo-TEM as in ref 40. A 4 μL droplet of aggregate aqueous solution was deposited on a Quantifoil (Quantifoil Micro Tools GmbH, Germany) holey carbon grid. The excess of liquid on the grid wasabsorbed with filter paper, and thegrid wasquench-frozen quickly in liquid ethane to forma thin vitreous ice film. Once placed in aGatan 626 cryo-holder cooled with liquid nitrogen, the samples were transferred to the microscope and observed at alowtemperature(-180°C).Observationswereperformed with aLaB6 JEOL JEM2100cryo-microscopeoperatingat 200 kV (IMPMC, Sorbonne Universite, CNRS UMR 7590, Paris, France) coupled with a 2k × 2k Gatan Ultrascan 1000 CCD camera. Images were taken with the JEOL low dose system (Minimum Dose System, MDS) to protect the thin ice film from any irradiation before imaging and to reduce the irradiation during the image capture.

UV-visible spectra were obtained in an Enspire multimode plate reader (PerkinElmer) using TPP 96-well plates in the 300-800 nm range.

Photothermal measurementswereperformed with a680 nm laser (Laser components SAS) illuminating a 0.5 mL tube containing0.1mL of sample.Thedistancebetween thesample and the laser was fixed to 4 cmto get a laser power of 6.3 W• cm -2 . The temperature of the sample was measured using a thermal camera (Flir SC7000) in real time (1 image/s) in the temperaturerangefrom20to 80°C.Datawereprocessed with the Altair software.

Theoretical Methods. All calculations were performed on R software (version 3.4.3).

To discriminate the different types of aggregates, the samples were represented using the PCA tool of the package FactoMineR and clusterized using k-means clustering.

TheSAR wasmodeled bylinear modelingusingthefunction lm of R. In the present case, the use of nonlinear terms was added to the lm function. All of the models were first calculated using the widest number of descriptors; then, to minimize the number of descriptors, a stepwise approach was performed on the model (function step, direction = "both", penalty fixed to log(n) for the linear models, 2log(n) for the nonlinear ones).

For the training dataset, the model quality was validated by assessment of residue normality (Figure S18), assessment of the equal distribution of residues through condition (Figure S19), and Y-randomization (R Y-random 2 = 0.10). Internal validation wasperformed byLeaveManyOut CrossValidation (LMOCV), usingthe package caret of R (function train). The model parameters were calculated on 50%of the data. Then, the accuracy of these parameters was evaluated by predicting R CV 2 on theremaining50%of thedataset.TheR CV 2 havebeen averaged on 1000 partitions of the data set. For the validation data set, the model prediction was validated by assessment of the equal distribution of residues through conditions (Figure S20).

The model accuracy on the validation data set was assessed with the function predict. The optimized descriptors to maximize the SAR were determined using the function optim. Asthefunction presentsasinglemaximum, thestarting point has no incidence. The domain of definition for this maximum research was established using either the physical limitation of the descriptor or the minimal and maximal values observed in the data set.

Figure 1 .

 1 Figure 1. Characterizationsof theGNP aggregatesobtained 5 min after mixingfor onesizeof GNPs(D =30 nm) and onetypeof polyelectrolyte (chitosan 300kDa) at varyingconcentration ratio:electrophoretic mobility(A) and diffusion relaxation timeDLS t1 (B) obtained through LDE and DLSmeasurements,respectively;fractal dimension(C) deduced fromSAXSmeasurements.Cryo-TEM (D) wasperformed for thesamplethat was the closest to the neutrality ratio (see arrow in part A). Optical extinction (E) and SAR (F) were measured, respectively, by UV-visible spectroscopy and photothermal measurements (λ =680 nm, 6.3 W•cm -2 ).

Figure 2 .

 2 Figure 2. Principal component analysis(PCA) and k-meansclusteringreveal theexistenceof threebehaviorsamongthesamples(A). Correlation between descriptors and the two first principal components (B). SAR across clusters (C).

Figure 3 .

 3 Figure 3. Fitted valuesvsobserved valuesfor theSAR for thetraining (orange squares) and validation data set (blue triangle). The red dotted line indicates the ideal case. R 2 is the correlation coefficient.

Figure 4 .

 4 Figure 4. Estimation of the parameter weight in the validation data set. Calculations were performed by calculating the parameter values on the training data set and then by an estimation of the determination coefficient on the validation data set.

Figure 5 .

 5 Figure 5. Validation of the predicted optimum of the model. Plot of the SAR values depending on the fractal dimension (A) and on the concentration ratio between polyelectrolyte and gold (B). The red dotted rectangles indicate the area around the optimized value for the fractal dimension and aggregate size according to the regression.

Table 1 .

 1 Summary of the Descriptors of GNP Aggregates

	descriptor	method	signification	unit
	DiamNPs		gold nanoparticle diameter nm
	MassPoly		mass average molar mass	kDa
			of the polyelectrolyte
	ChargePoly		number of charge per	e•mol•g -1
			molar mass
	RatioPolyGold		mass concentration ratio	no unit
			between polyelectrolyte
			and gold -log scale
	Time		time between mixing and	h
			measure
	EM	LDE	electrophoretic mobility	m 2 •V -1 •s -1
	DLS t1	DLS	first relaxation time measured by DLS -	μs
			log scale
	DLS A1 DLS t2	DLS DLS	proportion of objects associated to DLS t1 second relaxation time measured by DLS -	no unit μs
			log scale
	DLS A2	DLS	associated to DLS t2 proportion of objects	no unit
	FD	SAXS	fractal dimension	no unit
	UVvis m1	UV-visible spectroscopy	plasmon band position of the first	nm
	UVvis w1	UV-visible spectroscopy	band width of the first plasmon	nm
	UVvis A1	UV-visible spectroscopy	plasmon band intensity of the first	no unit
	UVvis m2	spectroscopy UV-visible	plasmon band position of the second	nm
	UVvis w2	UV-visible spectroscopy	band width of the first plasmon	nm
	UVvis A2	spectroscopy UV-visible	plasmon band intensity of the second	no unit
	SAR	PTT	heating property of the	kW•g -1
			aggregate -log scale

ACKNOWLEDGMENTS A.B. received a PhD fellowship from the doctoral school Physique en Ile de France (EDPIF). V.M.-A. received a postdoc fellowship from the Association pour le Recherche contre le Cancer (ARC, Aides Individuelles, postdoctorant, dossier 20150603405). The authors thank the ANR CarGold-16-CE09-026, ANR Coligomere-18-CE06-0006, and the European Union's Horizon 2020 research and innovation programunder Grant Agreement No. 801305 for funding. We are grateful to the Imagoseine core facility of the Institut Jacques Monod (UMR7592, CNRS/UniversitéParis Diderot, Paris, France), member of the France BioImaging (ANR-10-INBS-65904) for the electron microscopy observations, Thomas Bizien and the SWING beamline team (Synchrotron SOLEIL, St-Aubin, France) for the SAXS experiments, Jean Michel Guinier (Sorbonne Universite, Museúm National d'Histoire Naturelle, CNRS, IRD, Institut de Mineŕalogie, de Physique des Mateŕiaux et de Cosmochimie (IMPMC), Paris, France), Alba Nicolas Boluda and Plan-Sangnier (Universitéde Paris, Laboratoire Matiere et Systemes Complexes, UMR7057, Paris, France) for help with photothermal measurements, and Pierre Bost (C3BI, USR3756, CNRS/Institut Pasteur, Paris, France) and Cyprien Gay and Gaelle Charron (Universitéde Paris, Laboratoire Matiere et Systemes Complexes, UMR7057, Paris, France) for fruitful discussion.

AUTHOR INFORMATION