Louis Brézin

Thomas Risler

Jean-Franc ¸ois Joanny

Spontaneous flow created by active topological defects

Topological defects are at the root of the large-scale organization of liquid crystals. In two-dimensional active nematics, two classes of topological defects of charges ±1/2 are known to play a major role due to active stresses. Despite this importance, few analytical results have been obtained on the flow-field and active-stress patterns around active topological defects. Using the generic hydrodynamic theory of active systems, we investigate the flow and stress patterns around these topological defects in unbounded, two-dimensional active nematics. Under generic assumptions, we derive analytically the spontaneous velocity and stall force of self-advected defects in the presence of both shear and rotational viscosities. Applying our formalism to the dynamics of monolayers of elongated cells at confluence, we show that the non-conservation of cell number generically increases the self-advection velocity and could provide an explanation for their observed role in cellular extrusion and multilayering. We finally investigate numerically the influence of the Ericksen stress. Our work paves the way to a generic study of the role of topological defects in active nematics, and in particular in monolayers of elongated cells.

Topological defects are singularities in the topology of the order-parameter field in a phase of broken continuous symmetry [START_REF] Mermin | The topological theory of defects in ordered media[END_REF][START_REF] Kleman | Points, Lines, and Walls: In Liquid Crystals, Magnetic Systems, and Various Ordered Media[END_REF]. They were first observed in nematic liquid crystals as early as 1904 by Lehman [START_REF] Lehmann | [END_REF] and sparked theoretical interest later by Frank and Kleman for their role in the mechanical and optical properties of liquid crystals [START_REF] Kleman | Points, Lines, and Walls: In Liquid Crystals, Magnetic Systems, and Various Ordered Media[END_REF][START_REF] Kleman | Defects in liquid crystals[END_REF]. In two-dimensional nematic liquids, topological defects are point-like singularities of the director-field orientation. They are classified by the defect strength or topological charge S, which is equal to the number of full rotations of the liquidcrystal director on a closed path of total angle 2π around the defect core [START_REF] Oswald | Nematic and Cholesteric Liquid Crystals: Concepts and Physical Properties Illustrated by Experiments[END_REF][START_REF] Harth | Topological Point Defects of Liquid Crystals in Quasi-Two-Dimensional Geometries[END_REF]. Due to the nematic symmetry, this topological charge is either an integer or a half-integer. The regions that surround a topological defect present strong gradients of the director-field orientation, with high energetic costs. Because of their relatively lower energy, the defects that are mostly observed in passive nematic liquid crystals have topological charge S = ±1/2, which is the lowest possible charge in absolute value.

Largely driven by biophysical applications, much attention has been attracted more recently to active matter [START_REF] Jülicher | Active behavior of the Cytoskeleton[END_REF][START_REF] Marchetti | Hydrodynamics of soft active matter[END_REF]. a e-mail: jean-francois.joanny@college-de-france.fr b e-mail: thomas.risler@curie.fr In active nematics, numerical simulations [START_REF] Giomi | Defect annihilation and proliferation in active nematics[END_REF][START_REF] Thampi | Velocity Correlations in an Active Nematic[END_REF] as well as analytical analyses [START_REF] Pismen | Dynamics of defects in an active nematic layer[END_REF] have revealed a very rich dynamics, dominated by the creation, motion, and annihilation of defects. Topological defects in living systems have been first observed in 1968 by Elsdale in monolayers of fibroblast cells [START_REF] Elsdale | Parallel orientation of fibroblasts in vitro[END_REF] and, since then, the topological charge of the defects has helped to identify the nematic or polar nature of cellular tissues [START_REF] Kemkemer | Elastic properties of nematoid arrangements formed by amoeboid cells[END_REF][START_REF] Duclos | Topological defects in confined populations of spindle-shaped cells[END_REF][START_REF] Saw | Biological Tissues as Active Nematic Liquid Crystals[END_REF][START_REF] Blanch-Mercader | Integer topological defects of cell monolayers -mechanics and flows[END_REF][START_REF] Blanch-Mercader | Quantifying material properties of cell monolayers by analyzing integer topological defects[END_REF]. Recent works looked at the influence of topological defects on the onset of threedimensional morphogenesis [START_REF] Maroudas-Sacks | Topological defects in the nematic order of actin fibres as organization centres of Hydra morphogenesis[END_REF][START_REF] Metselaar | Topology and Morphology of Self-Deforming Active Shells[END_REF]. Some material properties can also be extracted from the observation of the orientation of the director field around topological defects [START_REF] Blanch-Mercader | Integer topological defects of cell monolayers -mechanics and flows[END_REF][START_REF] Blanch-Mercader | Quantifying material properties of cell monolayers by analyzing integer topological defects[END_REF].

The properties of topological defects in active nematic liquid crystals are very different from the properties of their equilibrium counterparts. The active stress in an active nematic is proportional to the nematic orientational tensor. As regions surrounding a topological defect have pronounced nematic orientation gradients, they correspond to regions of strong active force density, which induce local flows around the defect. The structure of this flow depends on the topological charge of the defect, but defects in general are associated with vortices of the flow field [START_REF] Giomi | Defect annihilation and proliferation in active nematics[END_REF][START_REF] Thampi | Velocity Correlations in an Active Nematic[END_REF]. If the director field around the defect breaks a left-right symmetry along a particular direction, and in the absence of any external force, the flow created by active stress gradients drags the defect along that direction with a finite velocity, limited by viscous dissipation. Defects of topological charge S = -1/2 have a three-fold symmetry. Consequently, they have no preferred direction of motion and passively diffuse. On the contrary, defects of topological charge S = +1/2 break a 'head-tail' symmetry. In active nematics, they generically move along their head-tail axis, with a direction dictated by their contractile (motion toward the tail) or extensile (motion toward the head) nature.

Experimental observations on nematic cell monolayers show that topological defects are preferential sites of cell extrusion [START_REF] Saw | Topological defects in epithelia govern cell death and extrusion[END_REF][START_REF] Kawaguchi | Topological defects control collective dynamics in neural progenitor cell cultures[END_REF] and multilayering [START_REF] Sarkar | Crisscross multilayering of cell sheets[END_REF]. Similar observations have been made in bacterial monolayers [START_REF] Copenhagen | Topological defects promote layer formation in Myxococcus xanthus colonies[END_REF]. In a recent work by T. Sarkar and colleagues [START_REF] Sarkar | Crisscross multilayering of cell sheets[END_REF] on the spontaneous organization of cell monolayers on a solid substrate, the authors measured the flow patterns around topological defects of charges ±1/2. Some of the +1/2 topological defects were found to be pinned at specific locations on the substrate, while others had a spontaneous motion as expected in active nematics.

In the current paper, we study theoretically the flows generated by active, two-dimensional topological defects of charges S = ±1/2. One key ingredient of the hydrodynamics of nematic liquid crystals is the backflow, the hydrodynamic motion induced by the rotation of the director field [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF]. Analytical results on topological defects have been obtained by neglecting this backflow [START_REF] Giomi | Defect dynamics in active nematics[END_REF][START_REF] Rønning | Flow around topological defects in active nematic films[END_REF], while numerical studies have considered the coupling between the flow and the director fields in passive [START_REF] Toth | Hydrodynamics of topological defects in nematic liquid crystals[END_REF] as well as active [START_REF] Giomi | Defect annihilation and proliferation in active nematics[END_REF][START_REF] Giomi | Defect dynamics in active nematics[END_REF] nematic defects. We provide here analytical solutions that account explicitly for these effects. We compute the force necessary to stall an active topological defect of charge S = +1/2 and investigate the link between the two-dimensional flow patterns and the trigger of multilayering. We limit our study to a single, isolated defect in an infinite plane. Thereby, we rely on analytical solutions of the active hydrodynamic equations that depend on the bulk properties of the cellular tissue only. Performing a perturbation analysis around the passive defect orientation allows us to obtain analytical solutions for weakly active systems. We also estimate the size of the spatial domain over which this hypothesis is valid.

The paper is organized as follows: section 1 introduces the hydrodynamic equations of active nematics by giving the constitutive equations and computing the active force density generated by a topological defect of charge S = ±1/2. We then study in section 2 the dynamics of the defect, assuming that the nematic orientation is unperturbed compared to that of a passive defect of the same charge. We focus on +1/2 defects, but computations for -1/2 defects are detailed in Appendix E. We compute the flow field and the selfadvection velocity of the defect in section 2.1. In section 2.2, we compute the force necessary to pin a topological defect to a specific location on the substrate. In section 2.3, we study the effect of inhomogeneous cell divisions and cell deaths or extrusions on the self-advection velocity of the defect, in the case where the rate of cell production is linearly coupled to the differential pressure field. In section 3, we address the role of the rotational viscosity, which induces feedback of the flow on the nematic orientation. We first ignore the elastic Ericksen stress and set the flow-alignment parameter to zero in section 3.1. This allows for the computation of the flow and the self-advection velocity to linear order in activity analytically. Within the same hypothesis, the crucial role of the backflow on the stall force is detailed in section 3.2. The first-order correction to the director orientation as well as the spatial domain of validity of our calculations are discussed in section 3.3. Finally, we investigate numerically the role of the Ericksen stress on the self-advection velocity in section 3.4.

Hydrodynamic equations

We consider an active nematic deep in the nematic phase where the order parameter has a constant modulus. To describe its dynamics, we make use of the hydrodynamic theory of active nematic gels [START_REF] Marchetti | Hydrodynamics of soft active matter[END_REF][START_REF] Kruse | Generic theory of active polar gels: A paradigm for cytoskeletal dynamics[END_REF][START_REF] Joanny | Hydrodynamic theory for multi-component active polar gels[END_REF]. This very general framework is based on symmetries and conservation laws [START_REF] Joanny | Active gels as a description of the actinmyosin cytoskeleton[END_REF] and has been shown to adequately describe monolayers of elongated cells that self-organize into a nematic phase [START_REF] Saw | Biological Tissues as Active Nematic Liquid Crystals[END_REF][START_REF] Duclos | Perfect nematic order in confined monolayers of spindle-shaped cells[END_REF][START_REF] Duclos | Spontaneous shear flow in confined cellular nematics[END_REF][START_REF] Alert | Physical Models of Collective Cell Migration[END_REF]. More specifically, we aim at describing the dynamics of nematic topological defects as observed experimentally in monolayers of myoblast cells [START_REF] Sarkar | Crisscross multilayering of cell sheets[END_REF]. Because we are interested in monolayers of cells at confluence, we consider the active nematic as a one-constituent dense phase, with a fixed cell density.

Within the active hydrodynamic framework, the monolayer is described by a velocity field v and a director field p. The two-dimensional Frank-Oseen free energy F associated to gradients of the director-field orientation reads, in the one-constant approximation where the splay and bend constants are equal [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF]:

F = dx dy K 2 Tr (∇ ⊗ p) • (∇ ⊗ p) T - 1 2 h 0 p 2 , (1)
where ⊗ denotes the tensorial product, and h 0 is a Lagrange multiplier to ensure that the director is a unit vector. The evolution of the velocity and director fields is described by two coupled vectorial equations. The first equation is the force-balance equation, which reads, at steady state:

∇ • σ σ σ = 0 . (2)
Inertia has been neglected since the Reynolds number for typical cellular systems is much smaller than one.

The total stress tensor is further split into a passive σ σ σ p and an active σ σ σ a contribution. The passive contribution is given by the hydrodynamic theory of passive nematics [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF][START_REF] Martin | Unified Hydrodynamic Theory for Crystals, Liquid Crystals, and Normal Fluids[END_REF][START_REF] Forster | Hydrodynamics of Liquid Crystals[END_REF] as:

σ σ σ p =2η ũ + ν 2 (h ⊗ p + p ⊗ h -(p • h) 1) + 1 2 (h ⊗ p -p ⊗ h) + σ σ σ E -P 1 . (3)
Here, η is the shear viscosity; ũ is the traceless part of the symmetric part of the velocity-gradient tensor u = (∇ ⊗ v + (∇ ⊗ v) T )/2; ν is the so-called flow-alignment parameter; h = -δF/δp is the orientational field; σ σ σ E is the traceless part of the Ericksen stress tensor σ σ σ E , which generalizes thermodynamic pressure for anisotropic systems [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF][START_REF] Chaikin | Principles of Condensed Matter Physics[END_REF]; and P is the pressure, which is a Lagrange multiplier that ensures the incompressibility condition ∇ • v = 0, with 1 the identity tensor. The elastic force density associated with the Ericksen stress is obtained from the Gibbs-Duhem relation [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF][START_REF] Joanny | Hydrodynamic theory for multi-component active polar gels[END_REF][START_REF] Jülicher | Hydrodynamic theory of active matter[END_REF]]

∇ • σ σ σ E = -(∇ ⊗ p) • h . (4)
The active contribution is linked to the existence of a local cell alignment along the direction of the director field p. It reads

σ σ σ a = -ζ∆µ p ⊗ p - 1 2 , (5)
where ζ∆µ is a scalar quantity that measures the activity of the system. It is the product of a difference in chemical potential ∆µ, which is positive, and an Onsager coefficient ζ, which can have either a positive or a negative value. A positive value of ζ corresponds to an extensile active stress, in which the cells push along their long axis. The stress ζ∆µ is negative in a contractile system where the cells pull along their long axis.

The second equation describes the evolution of the director field p:

Dp Dt = 1 γ h -ν u • p . (6)
Here, Dp/ Dt

= ∂ t p + (v • ∇)p + ω ω ω • p is the co-moving co- rotational derivative of the director, where ω ω ω = (∇ ⊗ v - (∇⊗v) T )/2
is the vorticity tensor, and γ is the rotational viscosity. The orientational field h = -δF/δp is associated to changes of the free energy with respect to the director p. It is parallel to p in a non-flowing steady state. It is convenient to introduce the components of the orientational field, parallel and perpendicular to p, h and h ⊥ , respectively. The component h ⊥ controls the orientation of the director p, whereas the component h controls the modulus of p [START_REF] Kruse | Generic theory of active polar gels: A paradigm for cytoskeletal dynamics[END_REF].

In the following, we define the orientation angle ϕ of the director p as p = cos ϕ e x + sin ϕ e y in the fixed reference frame {e x , e y }. For passive nematics, the equilibrium configurations correspond to minima of the Frank-Oseen free energy given in eq. ( 1), which leads to ∆ϕ = 0. This equation has the solution ϕ(r, θ) = S θ for a topological defect of charge S, using polar coordinates (r, θ) centered at the defect singularity. Defects with the lowest absolute topological charge S = ±1/2 are the most stable. They are shown in fig. 1. In free monolayers of nematic tissue, they are the only observed defects.

For a director field corresponding to a passive defect, using the polar coordinates (r, θ) in the Cartesian basis of vectors {e x , e y }, the active force density reads:

f a = ∇ • σ σ σ a = - ζ∆µ 2r e x , S = +1/2 -cos 2θ e x + sin 2θ e y , S = -1/2 (7)
As represented in fig. 1, the active force density for a +1/2 topological defect is along the axis of symmetry. This drives a self-advection of the defect 'from head to tail' in the contractile case and 'from tail to head' in the extensile case. For -1/2 defects, because of the three-fold symmetry of the director field, the net total force on a disk of radius R vanishes, such that the default, even tough active, can only diffuse.

Topological defects in the limit of vanishing rotational viscosity

Throughout this section, we consider the case of vanishing rotational viscosity. In this case, the molecular field vanishes to satisfy eq. ( 6), so that δF/δp = 0 and the director orientation satisfies the equilibrium condition for passive nematics ∆ϕ = 0, illustrated in fig. 1. The force density associated to the Ericksen stress as given by eq. ( 4) vanishes. The Ericksen stress can then be ignored and the passive stress as given by eq. ( 3) reduces to σ σ σ p = 2η ũ -P 1.

Flow field and self-advection velocity of a +1/2 defect

Under this assumption, we now compute the velocity field of an isolated +1/2 topological defect in an infinite domain. We write the equations of motion in the reference frame of the defect. Force balance reads

η∆v -∇P -ξ(v + v 0 ) - ζ∆µ 2r e x = 0 , ( 8
)
where ξ is the friction coefficient of the cell layer with the underlying substrate and v 0 the self-advection velocity of the defect with respect to the substrate. This bulk equation is supplemented by two boundary conditions on the velocity field: it vanishes at the core of the defect in the reference frame of the defect v(0, θ) = 0, and it vanishes at infinity in the reference frame of the substrate v(∞, θ) = -v 0 .

Considering the incompressibility condition ∇ • v = 0, eq. ( 8) is solved by introducing the stream function ψ(r, θ), defined by v = ∇ × (ψ e z ), where e z = e x × e y . We further adimensionalize the equations using the length scale L = η/ξ and the time scale τ = 2η/(|ζ|∆µ). In the following, all dimensionless quantities are denoted with a tilde. Taking the curl of eq. ( 8), the dimensionless stream function ψ = (τ/L 2 )ψ satisfies

∆ [∆ ψ -ψ] = -s sin θ r2 , (9)
where r = r/L, s = sign(ζ) and ∆ denotes the Laplace operator with respect to the reduced variable r. A +1/2 topological defect is symmetric with respect to the x-axis. Therefore, the radial component of the velocity field is an even function of the polar angle θ, and the stream function an odd function of that variable. Given that eq. ( 9) is linear in ψ with a forcing term in sin θ, only the Fourier mode n = 1 contributes non trivially to the solution, and the stream function can be written as ψ(r, θ) = ψ(r) sin θ. The dimensionless velocity field has then polar components ṽr (r, θ) = ṽr (r) cos θ and ṽθ (r, θ) = ṽθ (r) sin θ, which can be derived directly from the generic solution of eq. ( 9) for ψ, as detailed in Appendix B.1. Imposing that the solution is non-divergent at infinity and vanishing in r = 0, the solution for ṽr and ṽθ read:

ṽr (r) = - s r I 1 (r) +∞ r K 1 (u)(- π 4 u 2 + u) du + K 1 (r) r 0 I 1 (u)(- π 4 u 2 + u) du (10) ṽθ (r) =s I 0 (r) - I 1 (r) r +∞ r K 1 (u)(- π 4 u 2 + u) du -K 0 (r) + K 1 (r) r r 0 I 1 (u)(- π 4 u 2 + u) du , (11)
where I 1 (u) and K 1 (u) are the modified Bessel functions of the first and second kind, respectively [START_REF] Abramowitz | Handbook of Mathematical functions[END_REF] (see Appendix A). The self-advection velocity can then be computed. Taking the limits of eq. ( 10) at infinity leads to

v 0 = - π 4 × ζ∆µ 2 ξη e x , (12)
in physical units. This result is similar to that of ref. [START_REF] Rønning | Flow around topological defects in active nematic films[END_REF]. A representation of this velocity field as well as of its amplitude as a function of r is shown in fig. 2a,c. Corresponding expressions for -1/2 defects are given in Appendix E.1 and plotted in fig. 2b,d. The solution [START_REF] Thampi | Velocity Correlations in an Active Nematic[END_REF] depends only on intrinsic properties of the system and can therefore be directly compared with experimental measurements to determine intrinsic parameters in experimental systems. Because of hydrodynamic screening, the velocity field decays as 1/r at infinity in the substrate reference frame, as does the active-force density. The total active force diverges linearly with system size but is balanced by the total friction force. These divergences exist here because we consider a single defect in an unbounded space. For a pair of defects of opposite charges +1/2 and -1/2, for example, the activeforce density and the velocity field decrease as 1/r 2 . The pressure field is obtained from the force-balance eq. ( 8):

P(r, θ) = - ζ∆µ 2 cos θ + P 0 , (13)
where P 0 is a reference pressure. The active stress creates a differential pressure between the head and the tail of the defect. In the contractile case (ζ∆µ < 0), the pressure is higher at the tail of the defect, and it is the opposite in the extensile case.

2.2 Stall force of a +1/2 defect Sarkar et al. observed two classes of +1/2 topological defects in confluent monolayers of C2C12 mouse myoblasts [START_REF] Sarkar | Crisscross multilayering of cell sheets[END_REF]: motile topological defects with a flow characteristic of selfadvected defects, as presented in sec. 2.1, and stalled, nonmotile defects. These defects where not immobilized by jamming, as it can occur in other cellular systems [START_REF] Garcia | Physics of active jamming during collective cellular motion in a monolayer[END_REF], as significant flows were observed around them. In addition, a oneto-one correspondence was observed between non-motile defects and the initiation of a second layer of cells due to extrusion of cells from the first layer, a process referred to as multilayering. Immotile +1/2 defects in an active monolayer can only exist in the presence of an external force that balances the total active force, other than the friction drag that structurally vanishes for a vanishing cell-velocity field. We determine in this section the local external force applied on the defect core necessary to pin an active, +1/2 defect. We introduce a core region of the defect of finite size Vortices around +1/2 defects have a rotational direction that depends on the sign of ζ. The three-fold symmetry for -1/2 defects does not create self-advection contrary to +1/2 defects. c,d Dimensionless azimuthal (orange) and radial (blue) velocities as a function of the dimensionless distance to the core r in the reference frame of the defect. For +1/2 defects (c), this velocity vanishes at the core and tends to a finite limit at infinity. For -1/2 defects (d), this velocity vanishes both at the core and at infinity.

a of the order of the correlation length of the nematic order [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF], which we consider as a region where there is no nematic order. We therefore model the monolayer as divided into two regions: perfect nematic order with p 2 = 1 as in section 2.1 for r > a, and isotropic with p 2 = 0 for r < a. This assumption allows for an analytical treatment and is justified when looking at an isolated defect over large distances compared to the characteristic length L: since the total active force diverges with system size, the active force that would result from the weak nematic order inside the finite core region is negligible compared to the active force from the bulk region in a large system. Indeed, recent work show that, for isolated defects, the self-advection velocity given by eq. ( 12) is valid even when considering the effect of the weak nematic order inside the core [START_REF] Rønning | Flow around topological defects in active nematic films[END_REF].

Let us call f =f e x the external force applied to the core region necessary to stall the defect. We hypothesize further that, in the core region r < a, the external force-density has the same mathematical form as the active force in eq. ( 8) to allow for a similar treatment. Force balance then reads:

η∆v -∇P -ξv - ζ∆µ 2r e x = 0 r > a (14) η∆v -∇P -ξv - f 2πar e x = 0 r < a . ( 15
)
At the boundary r = a between these two regions, the global flow and stress fields are continuous. This boundary condition leads to an analytical treatment under the hypothesis that the radius a of the core region is much smaller than the hydrodynamic screening length L = η/ξ, as detailed in Appendix C.1. To leading order in a/L, we obtain

f = η ξ π 2 ζ∆µ log a ξ/η . (16)
A simple way to understand the result of eq. ( 16) is to consider the defect velocity in an active nematic submitted to a force localized at the core of the defect. As the hydrodynamic equations are linear, the velocity is the sum of two contributions: an active contribution v 0 proportional to the active stress ζ∆µ and a passive contribution proportional to the applied force f:

v defect = v 0 + χf , ( 17
)
where χ is the two-dimensional mobility of the defect. To linear order, this mobility is that of a defect in a passive nematic. A naive guess is given by the Saffman-Delbrück result for the mobility of a disk of radius a in a two-dimensional fluid of shear viscosity η: 1/χ ∼ η/ log(L/a), where the large-and short-scale cutoffs are the screening length L and the defect core size a, respectively. Using the active selfadvection velocity given by eq. ( 12) leads to eq. ( 16), up to a numerical prefactor.

Effect of cell division/extrusion on the self-advection of a +1/2 defect

We now look at the role of the non-conservation of the cell number on the active flow created by a topological defect. This non-conservation stems from cell-division and cell-death processes, as well as from cellular extrusion from the monolayer. We introduce an effective cell-proliferation rate k, which accounts for these three processes. With incompressible cells, the continuity equation then reads:

∇ • v = k . ( 18
)
There is evidence that the effective cell-proliferation rate is influenced by the cellular mechanical environment, and particularly by tissue pressure [START_REF] Basan | Homeostatic competition drives tumor growth and metastasis nucleation[END_REF][START_REF] Montel | Stress Clamp Experiments on Multicellular Tumor Spheroids[END_REF][START_REF] Delarue | Mechanical Control of Cell flow in Multicellular Spheroids[END_REF][START_REF] Delarue | Compressive Stress Inhibits Proliferation in Tumor Spheroids through a Volume Limitation[END_REF]. At a specific pressure called the homeostatic pressure P h , cell divisions and cell deaths and extrusions balance on average, such that k(P h ) = 0. Far from the defect, the tissue is in its homeostatic state P(r = ∞) = P h . In the vicinity of the homeostatic state, the pressure-dependent division rate reads to linear order:

k(P) = - 1 κ (P -P h ) , ( 19
)
where κ is a phenomenological coefficient. It is positive to ensure the monolayer stability and can be interpreted as an effective, long-term, bulk viscosity.

To solve the force-balance equation ( 8) in this context, we write the velocity field as a sum of a curl-free part and a divergence-free part, also known as the Helmholtz decomposition: v = ∇ × (ψ e z ) + ∇φ. Taking the curl and the divergence of eq. ( 8) leads to:

∆ [∆ ψ -ψ] = -s sin θ r2 ( 20
)
∆ ∆ φ - η η + κ φ = -s η η + κ cos θ r2 , (21)
using the same adimensional variables and functions as in Section 2.1. These equations are solved in Appendix D in a similar way as in sec. 2.1 to obtain the velocity field. The result for the self-advection velocity of the defect reads:

v 0 = - π 4 1 + η η + κ × ζ∆µ 2 ξη e x . (22)
The pressure-dependent division rate promotes motion of the defect in the same direction as what is dictated by the active force. We recover the incompressible limit of eq. ( 12) in the case of an infinite κ, that is when the cell number is conserved. In the other limit κ = 0, corresponding to an infinite response of the effective proliferation rate to pressure variations, the velocity of the defect is multiplied by two. In a contractile nematic, active forces associated with a +1/2 defect create a low-pressure environment at the head of the defect, which promotes cell division, and the opposite situation at the tail, as illustrated in fig. 3a. At large distances, the divergence of the velocity field vanishes, as illustrated in fig. 3c. This pattern of division and death/extrusion further induces motion form head to tail, similarly to the active force. For extensile systems, the situation is reversed. For a -1/2 defect, the pattern of effective cell divisions follows a three-fold symmetry. Corresponding expressions are given in Appendix E.2 and illustrated in fig. 3b,d. Interestingly, such a divergence profile has been observed in collective self-organizations of bacteria [START_REF] Copenhagen | Topological defects promote layer formation in Myxococcus xanthus colonies[END_REF].

Coupling between the velocity field and orientation gradients of the nematic director

All calculations in section 2 are performed in the limit of vanishing rotational viscosity γ. In this limit, the orienta- The director-field orientation is represented by solid, thin lines. In the contractile case represented here (ζ < 0), there is an asymmetric divergence profile for the +1/2 defect, with a positive divergence (orange) at the head and a negative one (blue) at the tail. For the -1/2 defect, there is a three-fold symmetry alternating between a positive (orange) and a negative (blue) divergence around the defect. c,d. Radial dependence of the dimensionless divergence for a +1/2 (c, ∇ • ṽ = D(r) cos θ) and a -1/2 (d, ∇ • ṽ-= D -(r) cos 3θ) defect. The amplitude of this quantity is maximum at the core of the defect and decays with the distance to the core.

tion of the director is obtained by minimizing the free energy eq. ( 1). Consequently, the velocity field was computed with the fixed director orientation ϕ = θ/2, corresponding to the orientation field of a passive defect at equilibrium. We now relax this hypothesis to look at the coupling between the flow field and the orientation field of the director. We therefore expect that activity influences the orientation of the director, further modifying the flow field with respect to the results of section 2. Since no global analytical solution can be found to this new system of equations, we compute the velocity field in perturbation in activity. We limit our investigation to first-order terms in activity. Since the flow field is already of order one in activity with the passive defect orientation ϕ = θ/2, the first-order modifications of the orientation contribute only to a second-order term in this velocity field. The first-order contribution of activity to the velocity field can therefore be computed using the passive defect orientation ϕ = θ/2. As a first approximation, we consider the case where the flow-alignment parameter ν vanishes, and we further neglect the Ericksen stress tensor σ σ σ E . The passive contribution to the stress tensor then reads

σ σ σ p = 2η ũ + 1 2 (h ⊗ p -p ⊗ h) -P 1 , (23)
and the molecular field h, given by eq. ( 6), reads at steady state:

h = γ [(v • ∇)p + ω ω ω • p] . (24)
3.1 Flow field and self-advection velocity

The active and passive contributions to the stress tensor as given by eqs. ( 5) and ( 23) lead to the force-balance condition:

η∆v + 1 2 ∇ × (h ⊥ e z ) -∇P -ξ(v + v 0 ) - ζ∆µ 2r e x = 0 . ( 25
)
The perpendicular component of the molecular field h ⊥ reads, using the equilibrium orientational field ϕ = θ/2:

h ⊥ = γ 2 v θ r -ω s , (26)
where ω s = (∇ × v) • e z is the vorticity. With this expression, the curl of the force balance equation ( 25) leads to:

∆ ∆ ψ - λ r ∂ r ψ -ψ = -s sin θ r2 , (27)
where we have adimensionalized the equation using the same procedure as in sec. 2.1, except with an effective viscosity η = η + γ/4, and where λ = γ/(4 η). More precisely, we have defined a new lengthscale L = η/ξ and a new timescale τ = 2 η/(|ζ|∆µ), defining the dimensionless variable r = r/ L and the dimensionless field ψ = (τ/ L2 )ψ. The term inside the brackets of the Laplace operator in eq. ( 27) is a transformed version of the modified Bessel equation [START_REF] Bowman | Introduction to Bessel Functions[END_REF].

Similarly to what was done in sec. 2.1 and as detailed in Appendix B.2, we determine the velocity field. The resulting self-advection velocity reads:

v 0 = -c(λ) × ζ∆µ 2 ξη e x , (28)
where

c(λ) = Γ [1 -β/2] Γ [1 + 1/(2β)] 2 √ λ Γ [3/2 -β/2] Γ [3/2 + 1/(2β)] (29)
with Γ the Euler gamma function [START_REF] Abramowitz | Handbook of Mathematical functions[END_REF] and

β = 1 + λ 2 /4+ λ/2.
We plot in fig. 4 the dimensionless amplitude ṽ0 of the self-advection velocity as a function of λ, both not considering and considering the Ericksen stress (see section 3.4). The rotational viscosity overall decreases this amplitude, which is expected since it corresponds to an additive source of dissipation. Fig. 4 Self-advection velocity of a +1/2 defect in units of ζ∆µ/(2 ξη) and as a function of λ. The dotted horizontal line shows the value π/4, corresponding to eq. ( 12) when λ = 0. Two cases are presented: without the Ericksen contribution to the stress (blue curve), as determined by eq. ( 28), and with the Ericksen contribution (orange curve), determined numerically (see section 3.4).

Stall force

Proceeding similarly as in section 2.2, the in-core force balance is unchanged as compared with eq. ( 15), since there is no nematic order within the core region, and the out-core force balance condition is given by eq. ( 25). Imposing continuity of the flow and stress fields at the boundary r = a between the in-core and out-core regions as detailed in Appendix C.2 gives, to leading order in a/ L:

f = - η ξ ζ∆µ g(λ) a ξ η β-1 , (30)
where

g(λ) =2 β-λ Γ [1 -β/2] Γ [1 + 1/(2β)] × 6π(2 -λ) [λ(4 -λ(1 + β)) + 24(β -λ/2 -1)] (48 -20λ 2 + λ 4 )Γ(1 + β -λ/2) . (31)
The force vanishes as a power law of the core size a with an exponent β -1, which varies between 0 and ( √ 5 -1)/2 as λ varies between 0 and 1. Compared to the previous case γ = 0, the stall force here is comparatively smaller in the limit of a small core size compared to the hydrodynamic length L, with a power-law dependence on the core size a rather than a logarithmic one.

First-order active correction to the director orientation

We can now investigate the effect of the flow on the orientation of the director, to first order in activity. With a non-vanishing rotational viscosity γ, the perpendicular component of the molecular field is given by eq. [START_REF] Rønning | Flow around topological defects in active nematic films[END_REF]. Given the definition of the molecular field and the expression of the free energy eq. ( 1), we obtain, in terms of the stream function ψ,

K∆ϕ = γ 2 ∆ψ - ∂ r ψ r . (32)
For reasons similar to those presented in sec. 2.1, the stream function remains of the form ψ(r, θ) = ψ(r) sin θ, and we obtain a similar form for the term proportional to the rotational viscosity γ in the expression of ϕ(r, θ):

ϕ(r, θ) = θ 2 + γ 2K {rv r (r) + 1 2 r r 0 v θ (u) u du - 1 r r 0 uv θ (u) du sin θ . (33)
Let us now investigate the behavior of the correction δϕ(r, θ) to the passive orientation ϕ 0 (r, θ) = θ/2 at large distances r L. Asymptotically, we get:

δϕ(r, θ) ≈ r L γ 4K r ln r r 0 v 0 , (34)
where r 0 is an arbitrary length that has been introduced based on dimensional analysis. Since this expression diverges with r, the solution given by eq. ( 33) is only valid up to a characteristic lengthscale L defect ∝ K/(γv 0 ). Using the scaling of eq. ( 28), we define

L defect = ξ(η + γ/4) γ/4 K |ζ∆µ| = L 2 a λ L , (35)
where L a = K/|ζ∆µ| is the active length, corresponding to the scale above which uniform bands of active nematics become unstable [START_REF] Voituriez | Spontaneous flow transition in active polar gels[END_REF][START_REF] Edwards | Spontaneous flow states in active nematics: A unified picture[END_REF] and active turbulence is expected [START_REF] Giomi | Geometry and Topology of Turbulence in Active Nematics[END_REF][START_REF] Doostmohammadi | Active nematics[END_REF]. This result gives an insight into the interplay between active and hydrodynamic stresses. Note that our description is only valid if L a L. In this limit, eq. ( 34) is valid at distances r from the core region such that L r L a . If this is not the case, the homogeneous nematic system does not have a stable steady state and the passive orientation ϕ 0 (r, θ) = θ/2 is spontaneously destroyed over length scales of order L.

Contribution of the Ericksen stress tensor to the self-advection velocity

Until now, we have only considered the coupling between the velocity field and the nematic order through the torques associated with the antisymmetric component of the stress tensor. In the expression of the stress tensor, we have however neglected the Ericksen stress tensor σ σ σ E as well as the flow-alignment coupling term proportional to the parameter ν. In perturbation theory around a homogeneous steady state, neglecting the Ericksen stress is in general justified by the fact that it is second order in the gradients of the director orientation, and therefore nonlinear in activity [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF]. Here, however, the Ericksen stress is relevant because the ground state upon which activity sets in has a non-homogeneous director orientation ϕ = θ/2. The force density resulting from the Ericksen stress is given by the Gibbs-Duhem relation [START_REF] Kleman | Defects in liquid crystals[END_REF], which here reads

∇ • σ σ σ E = - h ⊥ 2r e θ . (36)
Accounting for the Ericksen stress while keeping ν = 0, the curl of the force-balance condition leads to

∆ ∆ ψ - λ r ∂ r ψ -ψ + λ r ∂ r ∆ ψ - ∂ r ψ r = -s sin θ r2 , (37)
with the same dimensionless units as in eq. ( 27). Contrary to the previous equations for the stream function in this paper, this equation has no direct analytical solution. Using asymptotic expansions close to the core and at large distances, we can however determine the velocity of the defect using a shooting method. We plot in fig. 4 the resulting self-advection velocity of the defect as a function of λ (orange curve), in comparison to the analytical solution obtained in the absence of Ericksen stress (blue curve): there is a decrease in amplitude due to the Ericksen stress, which is expected given that it represents an additional source of rigidity.

Discussion

In this manuscript, we have studied the flows generated by ±1/2 topological defects of an active nematic system in two dimensions, which interacts with its underlying substrate via viscous drag. The active nematic system lies deep in the nematic phase, with an order parameter of constant, maximal amplitude. Our study is inspired by monolayers of elongated cells at confluence. In the case of a vanishing rotational viscosity, we derived analytical expressions for the flow field.

In the case of a +1/2 topological defect, we determined its self-advection velocity and stall force. We showed that the non-conservation of cell number leads to an increase in the self-advection velocity. We then investigated the effect of the rotational viscosity, which couples the velocity field to the nematic order. Taking the flow-alignment parameter to zero and in the absence of the Ericksen contribution to the stress, we derived analytical expressions of the self-advection velocity and stall force, making use of the passive, equilibrium configuration of the director around the defect. We finally studied the validity domain of the small departure from the equilibrium nematic-order configuration, as well as the influence of the Ericksen stress.

Our results show that the self-advection velocity is linear in activity for a +1/2 defect, as it should, as well as its corresponding stall force. The self-advection velocity of a +1/2 defect decreases in the presence of a non-vanishing rotational viscosity, which adds an extra source of dissipation. Accordingly, the stall force is also smaller, with a qualitatively different dependence on the size of the core region of the defect, where the nematic order vanishes. On the contrary, the non-conservation of the cell number increases the self-advection velocity, when coupled linearly to the departure from a homeostatic isotropic stress. Finally, accounting for the Ericksen stress lowers even further that velocity, with a similar qualitative interpretation.

Several assumptions used in our study can be questioned. A first restriction consisted in studying an isolated, topological defect in an infinite domain, as defects in an actual, free monolayers come into pairs to ensure a vanishing total topological charge. This assumption allowed for an analytical treatment that does not depend on finite-domain boundary conditions. Importantly, section 3.3 provided a characteristic system size, within which our calculations are valid, and beyond which active effects drive the nematic far from its passive, equilibrium configuration. At a qualitative level, the presence of other defects introduces another length scale in the problem, namely the characteristic distance d between defects. If this distance is such that d L, our results remain valid as the defects interact only very weakly. The interactions between defects becomes relevant in the limit where d L. In this case, one should consider the nematic layer as a gas of interacting defects [START_REF] Shankar | The low noise phase of a 2d active nematic[END_REF]. Within our approach, computing the flow created by a pair of +1/2 and -1/2 defects is the next logical step. A second restriction was to consider only partially the coupling between flow and orientation. In particular, for a non-vanishing rotational viscosity, we have limited our study to the case of a vanishing flow-alignment parameter. Beyond this approximation, no analytical solutions to the flow equations could be found. Finally, considering proliferation and extrusion is of prominent relevance for living systems, and we provided a first computation of this effect on the self-advection velocity of +1/2 defects.

The initial goal of this work was to explain the observation of Sarkar et al. in ref [START_REF] Sarkar | Crisscross multilayering of cell sheets[END_REF] of motionless defects that are preferential sites for multilayer formation. We computed the self-advection velocity of motile defects and the force necessary to stall them. However, the origin of the stalling and the mechanism by which another layer of cells forms at the defect is still elusive. Accounting for pressure-dependent proliferation shows that, on average, more extrusion is to be expected compared to the rest of the monolayer at the tail of +1/2, contractile defects. Recent studies by Vafa & Mahadevan [START_REF] Vafa | Active nematic defects and epithelial morphogenesis[END_REF] and by Hoffman et al. [START_REF] Hoffmann | Defectmediated morphogenesis[END_REF] show that, when considering a deformable surface, there is an out-of-plane force at topological defects that could be responsible for extrusion and multilayering. The boundary condition v c (r = 0) = 0 imposes that B c , B c 0 , and A c 0 vanish. The other integration constants are set by imposing the continuity of the velocity and stress fields at the boundary of the core r = a. Introducing P, the pressure normalized by |ζ|∆µ/2 and taking the divergence of the forcebalance eqs. ( 14) and ( 15), we obtain:

∆ P = s cos θ r2 r > ã (C.

Fig. 1

 1 Fig. 1 Passive orientation of the director field (solid lines) and resulting active force density (arrows) in the contractile case ζ∆µ < 0, for +1/2 (a) and -1/2 (b) topological defects. The origin of the polar angle θ is set by e x .

Fig. 2 a

 2 Fig.2a,b Flow structure around active +1/2 (a) and -1/2 (b) topological defects. Thin, solid black lines represent the orientation of the director field. Thick, arrowed blue lines are the directed streamlines of the velocity field in the reference frame of the substrate for a contractile nematic. Vortices around +1/2 defects have a rotational direction that depends on the sign of ζ. The three-fold symmetry for -1/2 defects does not create self-advection contrary to +1/2 defects. c,d Dimensionless azimuthal (orange) and radial (blue) velocities as a function of the dimensionless distance to the core r in the reference frame of the defect. For +1/2 defects (c), this velocity vanishes at the core and tends to a finite limit at infinity. For -1/2 defects (d), this velocity vanishes both at the core and at infinity.

Fig. 3 a

 3 Fig.3 a,b. Colormap of the dimensionless divergence of the velocity-or equivalently of the net division rate k-around a +1/2 (a) and a -1/2 (b) defect. The director-field orientation is represented by solid, thin lines. In the contractile case represented here (ζ < 0), there is an asymmetric divergence profile for the +1/2 defect, with a positive divergence (orange) at the head and a negative one (blue) at the tail. For the -1/2 defect, there is a three-fold symmetry alternating between a positive (orange) and a negative (blue) divergence around the defect. c,d. Radial dependence of the dimensionless divergence for a +1/2 (c, ∇ • ṽ = D(r) cos θ) and a -1/2 (d, ∇ • ṽ-= D -(r) cos 3θ) defect. The amplitude of this quantity is maximum at the core of the defect and decays with the distance to the core.

I 1 K 1 I 1

 111 (u)(A c u 2 + f u + B c ) du -K 0 (r) + K 1 (u)(A c u 2 + f u + B c ) du . (C.29)

30 )

 30 solved by:P(r, θ) = A * r + B * rs cos θ (C.32) Pc (r, θ) = A c * r + B c * r -f cos θ , (C.33)

Note that special care must be taken when matching the fields inside and outside the core, because normalizations have been done in these two instances using two different characteristic lengths, L = η/ξ and L = (η + γ/4)/ξ, respectively.

Acknowledgements We thank P. Silberzan and T. Sarkar from Institut Curie for sharing their experimental results with us and useful discussions. LB received a PhD fellowship from the doctoral school Physique en Ile-de-France (EDPIF) and support by the Collège de France foundation and Institut Curie. This work received support from the grants ANR-11-LABX-0038, ANR-10-IDEX-0001-02.

The modified Bessel functions I α , K α are the two general solutions of the following equation:

A related equation to the modified Bessel equation (A.1) is given by Bowman [START_REF] Bowman | Introduction to Bessel Functions[END_REF]:

Two independent solutions of eq. (A.2) are r a I d (br c ) and r a K d (br c ).

Appendix A.2: Derivatives

The first-order derivatives of the modified Bessel functions are given by:

I 0 (r) = I 1 (r) (A.5)

We also use the following asymptotic expansions:

where γ denotes here the Euler's constant. At infinity, we use This appendix is dedicated to the full computation of the velocity starting from eq. ( 9). As justified in the main text, the stream function reads ψ(r, θ) = ψ(r) sin θ. Integrating one Laplace operator gives:

where A and B are integration constants. The radial dependence of the stream function then satisfies:

The homogeneous solution ψ0 to this equation reads:

where A 0 and B 0 are integration constants. The Wronskian associated to eq. (B.12) reads

which leads to:

The velocity field ṽ = ṽr (r) cos θ e r + ṽθ (r) sin θ e θ is obtained from the derivatives of the stream function: ṽr (r) = ψ(r)/r and ṽθ (r) =d ψ(r)/ dr. We obtain

Imposing a finite velocity at the origin leads to B 0 = 0 and B = 0, as a finite velocity at infinity leads to A 0 = 0. We then obtain ṽr (r = 0) = -(A + sπ/4) and ṽθ (r = 0) = A + sπ/4 at the origin, and ṽr (r = ∞) = -A and ṽθ (r = ∞) = A at infinity. Imposing ṽr (0) = ṽθ (0) = 0 yields A = -sπ/4, and imposing ṽr (r = ∞) = -ṽ0 and ṽθ (r = ∞) = ṽ0 yields A = ṽ0 .

Coming back to the physical units, these lead to the selfadvection velocity given by eq. ( 12).

Appendix B.2: Velocity field with a finite rotational viscosity

We start the computation from eq. ( 27) with the dependence ψ(r, θ) = ψ(r) sin θ. Integrating one Laplace operator gives:

The left-hand side of this equation is of the form (A.2) with a = λ/2, b = 1, c = 1, and d = 1 + λ 2 /4. Following Appendix A, two independent homogeneous solutions to eq. (B.18) are rλ/2 I α (r), rλ/2 K α (r), with α = 1 + λ 2 /4. Following a similar procedure as in Appendix B.1, the dimensionless ve-locity field has the form:

The boundary conditions still set A 0 = 0, B 0 = 0, B = 0, and A = ṽ0 , similarly as in Appendix B.1. The asymptotic expansion of the velocity field close to the core reads:

Contrary to the case of a vanishing rotational viscosity, the velocity field in the defect reference frame vanishes close to the core for any finite value of ṽ0 , since the exponent λ/2 + α -1 is positive. To set the value of ṽ0 , we must consider the tangential stress σ rθ . Using the angular dependence of the stream function and force balance, its dimensionless version reads

Since λ/2 + α -2 is negative, ṽ0 must equal -sC λ 1 /C λ 2 for this tangential stress to remain finite at the origin. This leads to eq. ( 28) in physical units.

Appendix C: Stall force

Appendix C.1: Limit of vanishing rotational viscosity In the limit of vanishing rotational viscosity, the stall force is determined by the velocity and pressure solutions to eqs. ( 14) and ( 15), together with the incompressibility condition ∇ • v = 0. Quantities defined inside the core r < a bear the superscript 'c'. Using the dimensionless units of sec. 2, the stream function satisfies:

with ã = a/L and f = f /(πaζ∆µ) the normalized core radius and overall force applied to the core. The solution for the velocity field outside the core is given by eqs. (B. [START_REF] Blanch-Mercader | Integer topological defects of cell monolayers -mechanics and flows[END_REF]) and (B.17) with A 0 = 0 and A = ṽ0 as in Appendix B.1. The velocity field inside the core however reads: A * = 0 (C.36)

The normal and tangential components of the stress read σ rr = 2η∂ r v r -P and

We impose the continuity of the velocity and stress fields at the boundary of the core in r = a in the limit of a small core size ã 1, in which case we can make use of the following asymptotic expressions: ṽc r (r)

with C 0 1 and C 0 2 given respectively by eqs. (B.23) and (B.24) with λ = 0. Finally, we obtain the stall force by imposing ṽ0 = 0. The result is given by eq. ( 16), to leading order in ã = a/L = a ξ/η.

Appendix C.2: Finite rotational viscosity

The procedure to determine the stall force with a finite rotational viscosity resembles that presented in Appendix C.1. Since we assume no nematic order within the core region, the fields inside the core are unchanged as compared with Appendix C.1. The pressure field outside the core does not depend on the rotational viscosity and is unchanged. The quantity that changes is the velocity field outside the core region, now given by eqs. (B.21)-(B.24). Close to the core, for ā < r 1, the velocity components have the following asymptotic expressions:

with the notations of Appendix B.2. Using the expressions eqs. (C.32), (C.33), (C.38), (C.40), (C.42), and (C.43), continuity at the core boundary leads to the stalling force given by eq. ( 30) for ṽ0 = 0 1 .

Appendix D: Cell division/extrusion

In section 2.3, we introduce the Helmholtz decomposition v = ∇×(ψ e z )+∇φ. In dimensionless units, φ satisfies eq. ( 21), which is solved by:

where δ = (η + κ)/η. The velocity then reads:

where ψ is given by eq. (B.15) as in Appendix B.1. The force-balance condition (8) then imposes:

Imposing that the divergence of the velocity field-or equivalently the net division rate k-does not diverge at the core nor at infinity, we find that A 0 φ , A 0 , B 0 φ , B 0 , B, and B φ must all vanish. At the center of the defect r = 0, the components of velocity then read ṽr (r = 0) = -( ṽ0 + s(π/4)(1 + 1/δ)) and ṽθ (r = 0) = ṽ0 + s(π/4)(1 + 1/δ). Imposing ṽ(r = 0) = 0 in the reference frame of the default yields the self-advection velocity given by eq. ( 22).

The divergence of the velocity field is then given by:

It is represented in fig. 3.

Appendix E: -1/2 defects Appendix E.1: Velocity field This appendix is dedicated to the computation of the velocity field around a -1/2 defect, plotted on fig. 2. We add the superscript '-' to denote the quantities associated to a -1/2 defect. Force balance reads:

η∆v --∇P -ξv --ζ∆µ 2r (cos(3θ) e r + sin(3θ) e θ ) = 0 .

(E.50)

Here the defect is immotile, such that the velocity field satisfies:

The curl of eq. (E.50) gives:

Using a normalization by a characteristic length L = η/ξ and a characteristic time τ -= 2η/(3|ζ|∆µ), we introduce the dimensionless stream function ψ-= (τ -/L 2 )ψ -and spatial variable r = r/L. We then have:

where s = sign(ζ) and ψ -is of the form ψ-(r, θ) = ψ-(r) sin 3θ.

Integrating one Laplace operator in this equation leads to:

Solving this equation leads to the following dimensionless velocity field: We derive in this section the divergence of the velocity field represented on fig. 3. The derivation resembles that for +1/2 defects as presented in section 2.3 with a pressure-dependent division rate given by Eq. [START_REF] Metselaar | Topology and Morphology of Self-Deforming Active Shells[END_REF]. The velocity field is decomposed into v -= ∇ × (ψ -e z ) + ∇φ -. The divergence-free and curl-free parts of the velocity respectively satisfy:

∆ ∆ ψ--ψ-= s sin 3θ r2 (E.60)


Ben Huang

Jianling Yue

Benhui Fan

Yu Liu
email: yu_liu@csu.edu.cn

Xiaozhong Huang

Vertical carbon nanotubes arrays with controlled morphology on silicon carbide fibers for electromagnetic wave absorption

Keywords: CVD, SiC f -CNTs hybrids, electromagnetic wave absorption

Silicon carbide fiber-Carbon nanotubes arrays (SiC f -CNTs) hybrids with controllable morphologies were successfully synthesized by a floating catalyst chemical vapor deposition method (CVD). As the reaction duration prolonged from 2 min to 15 min, the length of CNTs was elongated from 2 μm to 50 μm, and moreover the content of CNTs in the hybrids increased from 1.5 wt.% to 3.7 wt.%. By filling this chopped SiC f -CNTs (20 wt.%) hybrids into paraffin, the composite with 3 mm in thickness presented the minimum reflection loss (RL min = -56.2 dB) and the maximum effective absorption wideband (EAB, RL≤-10 dB) of 4.3 GHz (8.5 -12.8 GHz). The dielectric performance of this hybrid synthesized by optimized CNT's length and growth density has presented high potential in the application of electromagnetic wave absorption.

Introduction

Reduction of electromagnetic wave (EMW) pollution has been required in the certain military or civilization applications such as airplane electronic systems and radar air traffic control communications [START_REF] Shahzad | Electromagnetic interference shielding with 2D transition metal carbides (MXenes)[END_REF][START_REF] Zhang | A Review on Metal-Organic Framework-Derived Porous Carbon-Based Novel Microwave Absorption Materials[END_REF]. With the advancement of nanotechnologies, the development of high-performance absorbing media with broad absorption bandwidth, low density, and thin thickness has aroused considerable attention from material science engineers [START_REF] Shahzad | Electromagnetic interference shielding with 2D transition metal carbides (MXenes)[END_REF]. Based on the theory of the interactions between EMW and material's microstructure, a wide range of materials with various structural parameters have been developed to reduce EMW pollution as large as possible. In the huge number of candidates, SiC fibers (SiC f ) have been attractive in this field because of their large aspect ratios, high thermal/chemical stability, and relatively low density. In particular, its high aspect ratio characteristic benefit SiC f to form interpenetrated conductive networks that accelerate the consumption of EMW energy through conduction losses [START_REF] Wang | Hierarchical Carbon Fiber@MXene@MoS2 Core-sheath Synergistic Microstructure for Tunable and Efficient Microwave Absorption[END_REF][START_REF] Hou | Implications from Broadband Microwave Absorption of Metal-Modified SiC Fiber Mats[END_REF]. However, as isolators, the low dielectric constant of SiC f partially limits the further improvement in the application of EMW absorption. Hence, many researchers have attempted to improve the dielectric properties of SiC f by various treatments. The introduction of fillers such as carbon materials [START_REF] Shi | Effects of thermal oxidation on microwave-absorbing and mechanical properties of SiC f /SiC composites with PyC interphase[END_REF][START_REF] Xu | High-Performance Electromagnetic Wave Absorbing CNT/SiCf Composites: Synthesis, Tuning, and Mechanism[END_REF], SiC nanowires [START_REF] Mi | Dielectric response and microwave absorption properties of SiCw/SiCf composites derived from carbon fiber[END_REF], metals [START_REF] Wei | Enhanced microwave-absorbing properties of FeCo magnetic film-functionalized silicon carbide fibers fabricated by a radio frequency magnetron method[END_REF], etc., on the surface of SiC f has been considered as an effective route to improve the dielectric permittivity that has been widely studied.

Carbon nanotube (CNT) is an effective filler that has been used in SiC f based composites for improving EMW absorption because of its specific features like high electrical conductivity, one-dimensional morphology, and the compatibility, etc. [START_REF] Chen | High-Performance Epoxy Nanocomposites Reinforced with Three-Dimensional Carbon Nanotube Sponge for Electromagnetic Interference Shielding[END_REF]. However, CNT shows a dominant mechanism of reflection of electromagnetic waves due to its high conductivity, which makes it difficult to achieve good wave absorption performance. Moreover, CNT is easy to agglomerate during the preparation, which makes it difficult to use directly as an electromagnetic wave absorbing material. To overcome those problems, using of fibers as a carrier to design the distribution pattern of CNT is an important means.

Chemical vapor deposition (CVD) has been shown to be an effective method for growing CNT. Among them, growing CNTs arrays on the SiC f substrate by floating catalyst chemical vapor deposition can effectively improve the CNT's incorporation and consequently makes the composites to achieve higher EMW absorbing [START_REF] He | Growth of carbon nanosheets on carbon nanotube arrays for the fabrication of three-dimensional micropatterned supercapacitors[END_REF]. The floating catalyst CVD is a major way to grow CNT array on various substrates [START_REF] Zhang | Tailoring the Morphology of Carbon Nanotube Arrays: From Spinnable Forests to Undulating Foams[END_REF][START_REF] Qin | Influencing factors and growth kinetics analysis of carbon nanotube growth on the surface of continuous fibers[END_REF][START_REF] Wang | Low temperature growth mechanisms of vertically aligned carbon nanofibers and nanotubes by radio frequency-plasma enhanced chemical vapor deposition[END_REF]. The control of synthesized CNT morphology on the substrates is not only related to the physicochemical properties and geometrical shapes of the substrates but also the reaction conditions.

In this work, we have successfully synthesized the growth of CNT array with controllable length and growth density on the SiC f substrates that allow us to systemically study the morphology effects of CNT on EMW absorption. Meanwhile, to explore the influences of sample size and the content of hybrids, the paraffin samples filled with chopped hybrids were prepared for the dielectric measurement in the microwave frequency range. The minimum reflection loss (RL min ) and effective absorption bandwidth (EAB, RL≤-10 dB) were calculated based on the tested dielectric permittivity. It was found that the performances of EMW absorption were highly related to the length of CNT in the hybrids. By optimization of CVD reaction conditions, the hybrids with suitable CNT's length and growth density can provide an effective approach to fabricate the composites with high EMW absorption performance.

Experiments and characterizations

The SiC f (Boxiang New Material Co.) was desized at 600 °C for 120 min in the air.

Then, the treated SiC f was immersed in 0.02 g/ml ferrocene/xylene solution for 60 min. The floating catalyst CVD method was used to deposit the CNT on the SiC f . The reaction temperature of the CVD method was 700 °C with an atmosphere mixed with Ar, H 2 , and C 2 H 2, while 0.05 g/ml ferrocene/xylene solution was used as catalyst precursor [START_REF] Liu | In-situ Growing Carbon Nanotubes on Nickel Modified Glass Fiber Reinforced Epoxy Composites for EMI Application[END_REF]. The studied reaction durations were 2, 5, 10, and 15 min, respectively.

The composition and micromorphology of the samples were characterized by thermogravimetric analysis (TGA), Raman spectrometer, X-ray diffraction (XRD), scanning electron microscopy (SEM), and transmission electron microscopy (TEM).

The electromagnetic parameters of the samples were tested by a vector network analyzer (Agilent 8720ET) in the frequency range of 2-18 GHz. The dielectric loss (

) was determined as the ratio between the imaginary part ( ) and the real part ( ) of relative dielectric permittivity [START_REF] Huang | Enhancement of the microwave absorption properties of PyC-SiCf/SiC composites by electrophoretic deposition of SiC nanowires on SiC fibers[END_REF]. Debye relaxation process is another important mechanism for absorbing EMW of dielectric loss materials. In spite of the different origins, various types of polarizations at microwave and millimeter-wave ranges can be described in a similar qualitative way. Considering the dielectric loss characteristics caused by the relaxation process, the relationship of and can be described as the following equation:

Results and discussion

(

Herein, and represent static dielectric constant and optic dielectric constant, respectively. According to Eq. ( 1), every single semicircle in the curve of -was called the Cole-Cole semicircle, corresponding to the Debye dipolar relaxation, and one semicircle is associated with one Debye dipolar relaxation process. As shown in The EMW absorption properties of the samples were evaluated by RL curves based on the measured electromagnetic parameters. According to the transmission line theory, RL at different frequencies (2-18 GHz) and thicknesses (1-5 mm) could be calculated by the following equations [START_REF] Zhu | Synthesis and electromagnetic wave absorption properties of peanut shell-like SiC fibers[END_REF][START_REF] Huang | Enhancement of the microwave absorption properties of PyC-SiCf/SiC composites by electrophoretic deposition of SiC nanowires on SiC fibers[END_REF]:

(2) (3)
Where and represented the input impedance and the impedance of free space;

and were the relative magnetic permeability and relative dielectric permittivity; and d represented the frequency and the corresponding thickness, and was the speed of light in the vacuum. Generally, the high-performance of EMW absorption depended on two factors, attenuation capability (α) and impedance matching (Z). On one hand, α was usually used to measure the attenuation capability of EMW for the absorber which is described as [START_REF] Dong | Construction of MnO nanoparticles anchored on SiC whiskers for superior electromagnetic wave absorption[END_REF]:

The variation of α value with frequency was shown in Figure 7 and On the other hand, besides α, the impedance matching (Z) also influenced the EMW absorption performance and it was described as the ratio between the input impedance ( ) and the impedance of the free space ( ), i.e., . When Z in was equal to or close to Z 0 , Z was equal to or close to 1, EMW would enter into the interior of the material to the greatest extent, rather than being reflected at the interface. In addition, the destructive interference theory accounted for the RL min would shift to low frequency with the increase of the absorber's thickness. The relationship of the corresponding thickness ( ) and frequency ( ) could be explained by the quarterwavelength equation [START_REF] Huang | Enhanced microwave absorption properties of carbon nanofibers functionalized by FeCo coatings[END_REF]:

(n=1, 3, 5 ...) (5)
where λ was the wavelength of EMW and was the measured frequency. The values of RL, , and Z were shown in Figure 8. The RL min was achieved at the frequency of the and the corresponding thickness of the when Z was equal to or close to 1. The closer the Z value to 1, the larger RL min could be achieved, and then the higher ability to absorb EMW. of the fibers makes it easy to form a conductive network, which also promotes the attenuation of EMW. Secondly, there are abundant interfaces in SiC f -CNT composites, which can effectively cause the scattering of EMW so that results in a consumption of energy. Finally, a suitable impedance match by adjusting the CNT content is favorable for the entrance of EMW into the SiC f -CNT composite, instead of being reflected directly on the surface of the composite. Therefore, the synergistic effects concluded by the mentioned three points assure the high-performance microwave absorption performance of SiC f -CNT composite.

Conclusion

In this work, SiC f -CNTs hybrids with controlled CNT array morphology were successfully synthesized by the floating catalyst CVD method. By controlling the CVD reaction conditions, the length and growth density of CNTs could be adjusted, which allowed us to study the morphological effects of CNTs on EMW absorption.

Besides, the EMW absorption performances were also associated with the content of the hybrids in the paraffin and the thickness of the sample. By comparison with the dielectric properties of different samples, the best absorption performance was obtained by the SiC f -CNTs hybrid with 15 minutes CVD reaction and 20% mass ratio

Figure 1

 1 Figure 1 (a) SEM image of desized SiC f ; SEM images of CNTs morphology on SiC f

Figure 1 Figure 1 (

 11 Figure1showed the SEM and TEM images of the desized SiC f and the SiC f -CNTs

Figure 2

 2 Figure 2 TGA analysis of samples: (a) room temperature to 800 o C; (b) 500 to 630 o C. TGA was conducted from room temperature to 800 o C in the air to evaluate CNT's

Figure 3

 3 Figure 3 XRD spectra for samples.

Figure 4

 4 Figure 4 Raman spectra for samples in 1000-1800 cm -1 .

Figure 4 presented

 4 Figure 4 presented Raman spectra of the samples. As presented by Raman spectra, the D band and G band appeared at 1330 cm -1 and 1598 cm -1 which were the typical

Figure 5

 5 Figure 5 Dielectric properties of the SiC f and SiC f -CNTs hybrids in paraffin with a

Figure 5 (

 5 Figure 5(d), there is a small tail of a straight line at the end of the -curves of

Figure 6

 6 Figure 6 Reflection loss curves for different thicknesses of 20 wt.% SiC f and SiC f -

Figure 7

 7 Figure 7Attenuation constant of the samples with a mass ratio of 20%.

  Figure S4, the α value increased with CVD growth duration, indicating that the content of CNTs should be the key factor for attenuation. Moreover, the same tendency of α with varied CVD growth durations was displayed in both Figure 7 and Figure S4. The large values of α would lead to impedance mismatch and consequently, result in degraded EMW absorption performance. Therefore, a reasonable value of α by tuning CVD duration could determine the EMW absorption performance of SiC f -CNTs.

Figure 8

 8 Figure 8 RL, t m , and Z values of different CNTs CVD growth times of 10 min with a filler loading of 20 wt.%.

Table 1

 1 EM wave absorption properties of CNTs on fibers.According to the above analysis, the excellent EMW performance and the wide EAB achieved by the SiC f -CNT composites mainly stem from the synergistic effects of strong intrinsic attenuation capability, interfacial scattering and the proper impedance matching. Firstly, the loss of the composites has increased as the CNTs content increases because of the high conductivity of CNT. Meanwhile, the high aspect ratio

	materials	Filler loading (wt.%)	RL min (dB)	EAB (GHz)	Thickness (mm)	Ref.
	CNTs/AHFS	40	-46.83	7.2	2.5	[20]
	CNTs-SiC f /SiC	-	-24	-	-	[21]
	CNTs/carbon fiber	20	-43.46	3.5	2.0	[22]
	CNTs/carbon fiber	-	-42	2.7	2.5	[23]
	CNTs/Co/carbon	10	-53.5	8.02	2.0	[24]
	fiber

Acknowledgment

This work is supported by the National Natural Science Foundation of China (Grant No. 52002403).

in paraffin: RL min of -56.2 dB (3 mm) with an EAB of 4.3 GHz. The excellent performances of SiC f -CNTs hybrids display high potential in EMW absorber media applications with multi-function.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.