Self-phoretic Brownian dynamics simulations

Abstract A realistic and effective model to simulate phoretic Brownian dynamics swimmers based on the general form of the thermophoretic force is here presented. The collective behavior of self-phoretic dimers is investigated with this model and compared with two simpler versions, allowing the understanding of the subtle interplay of steric interactions, propulsion, and phoretic effects. The phoretic Brownian dynamics method has control parameters which can be tuned to closely map the properties of experiments or simulations with explicit solvent, in particular those performed with multiparticle collision dynamics. The combination of the phoretic Brownian method and multiparticle collision dynamics is a powerful tool to precisely identify the importance of hydrodynamic interactions in systems of self-phoretic swimmers. Graphic Abstract


Introduction
Computer simulation of active matter systems is currently a topic of intense scientific debate [1][2][3][4][5]. Active matter considers systems with at least one component able to draw energy from their environment in order to self-propel. Activity is an inherent property of most biological systems and recently a topic of growing interest for the investigation of synthetic active systems, with practical applications in fields such as microfluidics or microsurgery [6,7]. In this line, phoresis is one of the main physical principles employed for the design of synthetic active matter. Phoresis refers to the drift that Brownian particles experience in the presence of a solvent with an intrinsic gradient, which becomes self-propulsion when the gradient is locally generated at the Brownian particle surface. Artificial microswimmers with a locomotion based on phoretic effects behave therefore as passive colloids unless activated via thermal [8][9][10][11][12], electric [13][14][15][16][17], chemical [18,19], or magnetic [20][21][22] gradients.
The collective behavior of chemically propelled Janus particles showed aggregation behavior [18,23,24], and light powered micro-robots were observed to form living crystals [25][26][27]. The appearance of clustering and comet-like swarming structures was predicted by Brownian thermophilic active colloids [28,29]. The system dimensionality [30,31] and the presence and shape of hydrodynamic interactions have shown to play a relevant role on the collective behavior of such thermophilic swimmers [32,33].
Janus-like phoretic particles have already been investigated by various simulation approaches, although a e-mail: s.roca@fz-juelich.de (corresponding author) b e-mail: m.ripoll@fz-juelich.de (corresponding author) not really compared with each other. Some of the approaches are purely Brownian, and self-phoretic propulsion is accounted simply by a constant impulse [28,34,35], or even a constant acceleration in systems that are supposed to increase their temperature on time [29]. In the absence of an explicit solvent, phoretic interactions between particles have been considered with an additional term, which might, or not, be coupled to the self-propulsion term. Thermal fluctuations are most frequently considered, and in a few cases also hydrodynamic interactions which are non-specific and typically only a far field approximation [36]. However, none of these methods completely accounts that self-phoretic Brownian swimmers propel with a well-defined Péclet number when isolated, while in the neighborhood of others, their velocity and interparticle interactions need to adjust to the actual distribution of heat sources. Different types of approaches consider the presence of an explicit solvent, such that phoretic effects arise in the presence of temperature or concentration gradients. This is the case of simulations performed with molecular dynamics [37], or dissipative particle dynamics [38], or with the mesoscopic simulation approach known as multiparticle collision dynamics (MPC) [11,39]. With these approaches, the details of self-propulsion, intercolloidal phoretic interactions, and hydrodynamic interactions are not directly imposed or tuned, but a consequence of the solvent-colloid interaction, colloid shape, and solvent intrinsic inhomogeneities. Therefore, in the studies of collective properties of phoretic active systems, the effects of steric, phoretic, or hydrodynamic interactions occur all simultaneously, such that the contribution of each of them is most frequently not possible to be identified. In this way, the design of strate-gies to make them distinguishable is timely and highly desirable.
Here, we propose a modification of the standard Brownian simulation method for Janus dimers, in which the effect of a single phoretic force results into the selfpropulsion of the dimer and interparticle interactions are included but without any hydrodynamic interactions. We refer to this method as phoretic Brownian dynamics (Ph-BD). Furthermore, the precise values of self-propulsion velocity, intensity of the interactions, and Péclet number can eventually be closely mapped to those of the MPC simulations to allow for a fair comparison of the results obtained with both methods. On the other hand, we discuss another two simpler Brownian dynamics types of approaches for dimers. One with only self-propulsion, and another with a constant selfpropulsion and phoretic interaction. An example study of dense systems of self-thermophilic dimers is here performed. The comparison of these three Brownian methods also provides interesting conclusions about the interplay of phoretic attraction/repulsion, alignment, and motility-induced instabilities.

Phoretic Brownian dynamics (Ph-BD)
Janus particles are characterized by having two different surface compositions. This is also the case of dumbbell-like structures in which each bead is made of a different material. For the sake of simplicity, we focus here mostly in the case of thermophoretic colloids, but the procedure is almost equivalent for other phoretic cases, in particular catalytic or diffusiophoretic ones. Thermophoretic dimers are made by one bead which is assumed to be at a higher temperature than the environment, mimicking a material with high heat conductivity which can be locally heated. The second bead is therefore exposed to a significant temperature gradient and responds to it depending on its intrinsic surface properties. In this way, dimers with a ther- Fig. 1 Sketches of the propulsion direction of self-phoretic asymmetric dimers and phoretic interaction between dimer pairs, which is: a attractive for the case of colloids drifting up gradient, and b repulsive in the opposite case. c Sketch of the implemented forces in the Ph-BD model mophilic (or chemotatic) behavior propel toward the hot bead (see Fig. 1a), while dimers with a thermophobic (or antichemotatic) behavior propel against the hot bead (see Fig. 1b). The same effect also controls the interaction between swimmers. Two thermophilic dimers swimming close to each other fell the temperature gradient produced not only by their own hot bead, but also by the hot bead of the neighboring dimer. This means that thermophilic dimers are attracted to neighboring dimers, while thermophobic dimers are mutually repelled from neighboring dimers, as depicted in Fig. 1a,b, which also exerts certain torques on the dimers. Therefore, in order to model phoretic active systems in a realistic manner, the effect of propulsion and interparticle interactions has to be included in a unified manner. Phoretic effects are related to the temperature gradients which vary locally.
Since the aim here is to describe the motion of colloids at low Reynolds numbers, we start by considering the overdamped Langevin equation [40,41] where F i (r) is the total sum of forces acting on each particle i, and ξ i (t) is a random force with zero-mean ξ i l (t) = 0, and delta-correlated Gaussian x, y, z and i, j = 1, . . . , 2N s the particles under simulation, with N s the number of simulated dimeric swimmers. The friction coefficient, μ i , is considered to fulfill the Stokes-Einstein relation, μ i = C f πηs i , with s i the radius of the particle i, and η the fluid viscosity.
The numerical factor C f varies depending on colloid boundary conditions, and it is typically C f = 6 for stick and C f = 4 for slip boundary conditions [42]. The algorithm here used to integrate the motion equations was stochastic Euler. In general, the Euler algorithm has to be carefully considered due to its low precision and problems in most isothermal simulations [43,44]. Nonetheless, the precision of this algorithm is sufficient for the study here performed, and other algorithms can be easily employed for more extensive investigations.
The details of the interactions are then provided by the forces, which distinguish two types of particles: with i, j = 1, . . . , N s for the two beads of each dimer. The non-heated or phoretic bead is the one where the temperature gradient has a drift effect, which in this Ph-BD approach is considered in an effective manner by including a thermophoretic force F T (see Fig. 1c). Meanwhile, the hot bead is considered to be at a higher but constant temperature, such that does not feel any thermophoretic force.
Pairwise interactions are considered, first the two beads forming each dimer are linked by a strong harmonic force F H , obtained from the potential with the interparticle distance r ij = r i − r j , the harmonic constant κ H = 10 4 , used to strongly fix the beads equilibrium distance b as the sum of the beads' radii, b = s p + s h , with s p the radius of the phoretic bead and s h the radius of the hot bead. The relative dimensions of the dimer beads are defined by the radii aspect ratio, γ = s p /s h . Steric effects are accounted by the force F EV,j , by which all non-linked beads interact with excluded volume interactions taken into account with the potential where n determines the potential softness, ε = k B T relates to standard energy units, with k B the Boltzmann constant and T the average temperature, here both fixed to unity, defining the system units. The extra term on the right-hand side of the equation determines the repulsive character of the potential together with the cutoff radius, r c = 2 1/n σ, and we use here n = 24. The distance σ simply relates to the sum of the radius of the two interacting beads. The thermophoretic force exerted on the phoretic bead can be calculated as where α T is the bead thermodiffusion coefficient and ∇ ri T the gradient of temperature at the bead location. Note that α T is a material property, which can be arbitrarily modified or chosen to match a value determined by experiments or by simulations with explicit solvent, as will be shown later.
The corresponding Laplace equations need to be solved to obtain a good estimation of the temperature gradient. We consider here three important and reasonable simplifications, such that the Laplace equation can be analytically solved: i) Each hot bead center acts as point-like heat source with temperature T h at the bead's surface; ii) at a distance far enough the fluid reaches the average fluid temperature T , taken as the reference unit, i.e., T (r → ∞) = T ; iii) the effect of neighboring sources is considered to be additive. For each point-like source, all angular terms vanish due to the symmetry of the system such that ∇ 2 T (r) = 0, and the temperature at the center of each phoretic particle is then given by where r ij = |r i − r j | and r j are the hot bead's center positions. The expression in Eq. (7) corresponds to the gradient at the bead center. A more accurate estimation is to consider an effective value of the temperature gradient that considers the variation over the bead surface, for which Eq. (7) is integrated along the phoretic bead's diameter. For a phoretic bead of radius s p , placed at r i , and with a hot bead placed at distance r j , the integral limits are r ij −s p and r ij +s p , such that the temperature gradient can be approximated by Note that for an isolated swimmer the gradient is determined just by the linked hot bead, such that r ij = s p + s h is the only contributing term. For denser systems, the gradient takes into account the all neighboring hot heads, such that in the center of highly compact configurations the gradients eventually vanish and therefore also the thermophoretic force. The dimer velocity v s and the rotational diffusion D r are therefore not direct inputs of the model, but indirectly determined from other input values, mainly α T , ∇T , s p , and γ. The value of the module of v s is given by where both the dimer friction μ = C f πη(s h + s p ) and the temperature gradient ∇T = (T h − T )/(s h + 2s p ) depend on the hot and phoretic bead sizes. The axis direction of the swimmer n aligns with the direction of the temperature gradient direction considering also the sign of α T , which determines the direction of v s . The self-propulsion velocity can also be obtained from the simulations as v s = v · n. The rotational diffusion D r depends mostly on the particle size and aspect ratio γ and can be obtained by characterizing the longtime behavior of the mean-squared angular displacement, Δe 2 = (e(t) − e(t )) 2 , in simulations with equilibrium conditions; this is with T h = T . The resulting Péclet number can then be defined as Pe = v s /(D r s p ).

Other active Brownian dimer models
The method proposed in this manuscript, Ph-BD, differs from other approaches employed in the literature in the way that phoretic self-propulsion and interparticle phoretic interactions are coupled to each other. In order to better understand the relevance of this coupling, we propose two alternative methods. Self-propelled spherical colloids have been extensively investigated with the so-called active Brownian particle (ABP) model [34,45], which simply assumes a constant propulsion velocity in the particle main axis. The physical origin of the propulsion is not specified, such that it could be phoretic but also any type of biological specificity. We adapt this idea to the dimeric case by considering N s swimmers with two bounded monomers each, where the hot bead just follows Eq. (3), and the phoretic bead where the friction is that of the dimeric structure μ = C f πη(s h + s p ) and n is the orientation vector of the dimer. We here call this method the active Brownian multimer model (ABM). With this approach, there is no additional interparticle interactions, such that all apparent repulsions or attractions are consequence of the propulsion and/or steric interactions. The second approach includes also the effect of the phoretic interaction with a force as given in Eq. (6), but considering only the heat sources of neighboring hot beads, where the temperature gradient can be calculated with Eq. (7) or Eq. (8). We refer to this method the active Brownian multimers with phoresis model (ABM+ph). With this approach, the phoretic interdimer attraction (or repulsion) is in principle decoupled from the dimer propulsion since there are two different parameters control, i.e., v s and α T . There are approaches in which these, or very strongly related parameters, are independently varied [34,35], which cannot really correspond to a phoretic model since both self-propulsion and interparticle phoresis are simultaneously originated. Besides the fact that v s and α T should be related by Eq. (9) for thermophoresis, or an equivalent one for other phoretic phenomena, there is another relevant difference between Ph-BD and ABM+ph which is that in ABM+ph, the velocity of the particles is fixed, namely it does not depend on the position of the neighboring particles, while for Ph-BD both the velocity of the particle and the interparticle interactions are damped when various other swimmers are in the neighborhood, as accounted in the temperature gradient calculation. In order to better understand this effect, we focus here in the case that v s and α T are linked by Eq. (9).

Hydrodynamic self-phoretic model
The methods introduced until now consider steric, stochastic, and phoretic interactions, which means that hydrodynamic interactions (HI) have been disregarded.
Although in some cases this can be clearly justified, the effect of HI is frequently not known. In order to provide a tool that allows for a fair comparison, we consider now the method known as multiparticle collision dynamics (MPC) [11,39]. Multiparticle collision dynamics is here used to simulate the explicit solvent particles and their interactions [46,47], while molecular dynamics (MD) is employed for colloid-colloid and colloid-solvent interactions. This hybrid MPC-MD approach has already extensively proved to include both hydrodynamics and phoretic effects [48][49][50][51].
MPC method for the solvent The MPC method considers the solvent composed of N point particles of mass m performing alternate streaming and collision steps. During the streaming step, fluid particles translate ballistically for a certain time, h, the collision time, this is r k (t + h) = r k (t) + hv k (t). In the collision step, the particles are binned into cubic cells of side a, with a grid shift applied to the binning in order to restore Galilean invariance [52]. Interparticle interactions are treated within each of these cells, in which particles interchange linear momentum with all other particles in the same cell. Here, we employ the stochastic rotational dynamics collision rule, in which the momentum interchange is made rotating by an angle α the relative velocities to the center of mass around a random axis on the cell,  [53][54][55][56]. The comparison with specific solvents can be done via dimensionless numbers, mainly the Schmidt number, Sc = ν/D = 13, and the Prandtl number, Pr = ν/κ T = 5.3. While Sc is smaller than the value for water, Pr is quite close to it. These two values ensure that momentum transfer is faster than that of mass, providing an efficient way to include hydrodynamic interactions, and that the stability of local temperature gradients is also ensured.

Molecular dynamics
Fluid-colloid interactions are considered using molecular dynamics, with the equations of motion being integrated using the velocity Verlet algorithm [43,44,57]. The thermophoretic nature of the colloids is determined by the choice of the fluidcolloid interactions, for which we used a displaced Mielike potential, This potential is very similar to that in Eq. (5) with the introduction of the Δ and the C parameters. The bead size is now determined by s ≡ σ + Δ, where Δ can be understood as the size of a core with hard-sphere interactions and σ the size of an additional layer with repulsive potential interactions. In this work, we use σ = Δ and s p = 6 for the size of the phoretic bead. The extra term on the right-hand side of Eq. (12) is C = for repulsive interactions, which have proved to account for thermophilic colloidal behavior, and C = 0 for attractive interactions for thermophobic [58,59]. For these interactions, n = 3 is chosen to obtain a soft repulsive potential for the phoretic (philic) bead, whereas n = 24 is chosen for the heated particle and also for the attractive potential (phobic). The cutoff radius of the interactions is r c = 1.26σ + Δ for the repulsive potential and r c = 1.1σ + Δ for the attractive. Harmonic and excluded volume interactions are considered similarly as for the Ph-BD case with Eq. (4) and Eq. (5). In order to mimic the heating produced by laser illumination of partially gold-coated colloids [60], we have rescaled the temperature of the fluid within a small shell (of 0.08s h ) around the heated bead to T h > T , while cooling the average temperature of the whole system to T = 1, by means of a simple velocity rescale [55,60]. Unless otherwise specified, we use T h = 1.5. All colloid-colloid interactions have been implemented via Eq. (5); this is Eq. (12), with Δ = 0, σ = s and n = 24, with the interactions being cut at r c = 2 1/24 σ.
This method has been implemented on LAMMPS [61], where we have modified the "srd" package routine [62] to include the colloid-solvent potential interactions. The MD time step has been chosen as Δt = 0.01h, similar as in the Brownian simulations, and the mass M of the colloidal beads is chosen to make the colloids neutrally buoyant.

Parameters for the comparison MPC vs. Ph-BD
In order to perform a fair comparison of the methods with and without HI, we are interested in having systems as similar as possible. Some values are input parameters in the Brownian dynamics simulations and therefore very easy to match, such as the average temperature k B T = 1, or the fluid viscosity η = νρ = 7.9.
The numerical factor C f for the friction coefficient is fixed as C f = 3 in order to match the employed MPC-SRD algorithm without angular momentum conservation and slip boundary conditions [63,64]. Other parameters are not direct input and need to be more carefully considered. For a proper comparison, it is of importance that parameters chosen for the two simulations models result in matching self-propulsion velocity and the Péclet number of diluted swimmer dimers systems. For this, we need to characterize the simulated thermophoretic coefficient α T and rotational diffusion D r .
The thermophoretic coefficient α T of a spherical bead could in principle be determined in full hydrodynamic simulations with an external temperature gradient [58]. This would be, however, a too rough estimation, first because the constant and position-dependent gradients are different, and second because it is known that the proximity of the hot bead screens part of the phoretic interactions of the surrounding solvent and the colloid surface. A more adequate estimation can be done by measuring the self-propelled velocity of a single swimmer with Eq. (9) and relating it then to the thermophoretic coefficient α T . Figure 2 shows simulation Results for dimers with phoretic bead sp = 6. Circles (in blue) correspond to thermophobic dimers; triangles (in red) correspond to thermophilic dimers. Full symbols correspond to asymmetric dimers (γ = 3); empty symbols to symmetric dimers (γ = 1). Lines relate to linear fits to Eq. (9) for small gradients results for four types of dimeric swimmers, corresponding to thermophobic and thermophilic character, and to the symmetric (γ = 1) and asymmetric (γ = 3) geometries. Velocities are calculated as an average of 20 independent simulations, and the error bars are of the order of the symbol size. Simulations are performed at various temperature gradients, which are achieved by changing the temperature of the hot bead T h . The increase in the velocity is clearly linear for moderate gradients, which allows us to determine the value of the thermophoretic coefficient for all the investigated cases as shown in Table 1. For the largest temperature gradients, the velocities deviate from the linear behavior. This deviation from the Fourier linear behavior can be expected and here can also be related to the limit of the method for these temperature gradients. Note that the negative sign of α T is well established by convention and it refers the motion of the swimmer toward the heat source. The sign of this coefficient naturally induces the interdimer phoretic attraction for thermophilic dimers, and phoretic repulsion for thermophobic ones, as shown in Fig. 1a, b, such that no further assumption has to be made in this regard.
To perform simulations with Ph-BD and MPC of dimers with the same bead sizes, the same solvent input parameters, and the same α T seems then a good strategy to have a fair comparison between methods, only the Péclet number is left to be discussed. The values of D r for the specified values result to be close to 60% larger in the Ph-BD simulations than in those performed with MPC, for all the four investigated cases (see Fig. 3a for the case s bd p /s hi p = 1). This is because the rotational diffusion is not a parameter fixed in any of the two methods but a consequence of all the other parameters, such as friction, particle size, thermophoretic coefficient, or fluid particle interactions which are different in both methods. In order to modify the rotational diffusion coefficient without affecting Table 1 Thermophoretic coefficient αT of single thermophilic and thermophobic self-propelled dimers, with different aspect ratios γ = s h /sp, as obtained from MPC simulations. Values are obtained as a fit to the data in Fig. 2 Table 1. The dashed lines at unity indicate perfect agreement between Ph-BD and MPC simulations. Thick vertical gray line corresponds to the case with optimal agreement for both vs and Pe, which occurs for s bd p = 8 other characteristic values, it is possible to vary the overall bead sizes. Further simulations in equilibrium with dimers with different s p (and different s h to preserve γ) are performed, and the measured values of D r are shown in Fig. 3a. As expected, the results show a decay of D r for growing particle sizes, which is produced just by the thermal noise in the position update of the two linked dimer monomers. Using a different monomer size allows then the tuning of the rotational diffusion, but, for simulations out of equilibrium with a fixed value of T h , also modifies the temperature gradient at the phoretic monomer surface and therefore also the resulting self-propelled velocity. The solution to keep the same v s value is then to modify T h to keep the gradient constant when changing the monomer size. As a check of this principle, we perform simulations modifying both s p and T h for a given gradient and then measure the self-propelled velocity. The results are shown in Fig. 3b in comparison with those of the self-propelled velocity of the hydrodynamic simulations used here as an input, and the agreement is very good within the error of the measurements in all cases. Note that in order to keep the same v s value, to modify T h is equivalent to modify α T since it is the product of both which determines F T and v s , as shown in Eq. (6) and Eq. (9), respectively. The resulting Péclet number shown in Fig. 3c results in a growing trend with particle size, which is related to the variation in the rotational diffusion. From Fig. 3, it is also clear that the optimal value is given by Brownian simulations with s p = 8, such that all presented Brownian simulations are from now on carried with this value.

Comparative study for collective dynamics
In order to perform a comparative study of the Brownian methods, simulations of dimeric thermophilic swimmers are performed first with the three Brownian methods previously discussed. Ensembles of 200 dimers both asymmetric, γ = 3, and symmetric, γ = 1, have been studied for a quasi-2d confinement case. In principle, this refers to 3d slides of liquid in which the swimmers move on a plane, which for the Brownian dynamics simulations means that the motion occurs in two dimensions.
The configuration used to initialize the simulation has the dimers center of mass placed on a square lattice covering almost the whole simulation box, with a randomly chosen direction of the dimer axis. Initial order disappears very quickly in all cases. All simulations run for a time t ∼ 300τ b , with τ b the ballistic time of a swimmer, defined as τ b = s p /v s , and representative snapshots of the latest configurations are shown in Fig. 4. Asymmetric dimers propel forming small clusters; some of these clusters dissolve due to collisions with isolated dimers or pairs of dimers, and some other coalesce with other small cluster, forming larger and more stable clusters, as can be seen in Fig. 4a. Symmetric dimers show initially similar dynamics, although interestingly the large clusters do not become stable and also end up dissolving in this case, as can be seen due to the small-sized clusters in Fig. 4b. Snapshots in Fig. 4 correspond to simulations performed with the Ph-BD method; qualitative roughly similar results are also obtained with ABM and ABM+ph methods. In order to more precisely understand the involved mechanisms and the difference between the methods, the quantification of a dynamic quantity is employed.
We introduce here the calculation of bounding time τ c for both asymmetric and symmetric dimers at two  The bounding times of asymmetric dimers in Fig. 5a show to form stable clusters at both densities in simulations with phoretic attraction; this is with Ph-BD and ABM+ph simulations. Simulations without phoretic interparticle attraction show to saturate to a constant value, which curiously is the same for both simulated densities. This means that for these asymmetric dimers at these densities, self-propulsion is not enough to stabilize the clusters, and the consideration of the corresponding phoretic attraction stabilizes the clusters. The bounding times are slightly smaller for ABM+ph with respect to Ph-BD. Although the difference is not large, it indicates that diminishing the propulsion velocity of the dimers inside the cluster slightly increases its cohesion. The bounding times of symmetric dimers in Fig. 5b show that all three Brownian methods form just unstable clusters; several interesting conclusions can still be drawn from these results. The first one relates to the effect of density, which shows to increase the bounding time in all cases, although with different intensity. Density increases the probability of encounters which has two opposite effects since it enhances both clustering formation and its dissolution. For the symmetric dimers without phoresis (ABM), the effect is small but clear at not too small times. This is in contrast to the asymmetric case for which the difference is much smaller and even seem to disappear for large averaging times. We relate this difference to the particleinduced alignment when two dimers collide, effect that is much larger for symmetric dimers. Another clear conclusion is that Ph-BD enhances stability with respect to ABM+ph and that this effect increases with density. This is again related to the fact that the smaller propulsion velocity of the dimers inside a cluster for the Ph-BD cases increases their stability. The effect is though not straightforward to predict since it is not shown to be much larger for asymmetric dimer swimmers where the formed clusters are larger, than in the case of symmetric dimers where only small clusters would be affected. Curious is also the difference between ABM and ABM+ph for the symmetric dimers. The collective simulations here presented only analyze the thermophilic case, such the inclusion of phoresis includes an interparticle attraction, expected to translate into larger clustering affinity. This is indeed the case for the asymmetric dimers, but not for the case of symmetric dimers. Phoretic attraction combined with the non-adjusting self-propulsion velocity seems to induce additional alignment of the symmetric dimers, such that they became more prompted to swim away from the small nucleated clusters, producing this somehow counterintuitive effect. In other words, the fact that ABM dimers do not significantly change their orientation when colliding with others makes that in some cases they get stuck in configurations longer than in the presence of attraction, providing such structures with additional stability. For the φ = 0.3 symmetric density case, the ABM simulations have almost the same stability properties as those with Ph-BD, while for φ = 0.2, ABM simulations are even more stable than those with Ph-BD, where both self-propelled velocity and attraction diminish in the neighborhood of other dimers.
Simulations in the collective regime with the hydrodynamics phoretic model MPC are also performed in the previously discussed regimes. The most relevant conclusion is the occurrence of qualitative differences with the systems here presented, which, due to the fair comparison that these methods provide, can be related just to hydrodynamics. Detailed understanding of such results requires a detailed discussion of the shape of the hydrodynamic fields which will be presented elsewhere.

Summary and discussion
The Ph-BD method to perform Brownian dynamics simulations with a realistic inclusion of the phoretic self-propulsion and interactions is here presented. The main idea is that simply considering the well-known dependence of the phoretic force with the applied gradient properly couples the self-propelled velocity and the interaction between two or more particles. The Ph-BD method is here compared with other two simpler versions of BD simulations, which draws interesting conclusions illustrating the very subtle interplay of selfpropulsion, phoretic-induced attraction, repulsion, or orientation. Depending on the particle geometry, properties, and overall densities, these effects show that they can act together, or against each other. We also show in detail how the Ph-BD method can be adjusted to map the properties of experimental systems or sim-ulations with explicit solvent, as illustrated here for MPC simulations. The combination of MPC and Ph-BD simulations offers therefore the possibility to compare simulated systems which only differ in the inclusion of solvent-mediated interactions, which is a very powerful tool to understand the effect of hydrodynamic interactions. These methods are here used to investigate the properties of thermophoretic dimers, but it can be almost trivially extended to other phoretic effects such as diffusiophoresis, and also to other multimeric structures such as trimers or other oligomeric swimmers. Preliminary analysis indicates that it can also be extended to Janus spherical particles. The presented results for the collective dynamics of thermophoretic swimmers also indicate that these are the basis of synthetic active materials with various perspectives for applied materials.