J Perdigon
email: jeremy.perdigon@oca.eu

G Niccolini

M Faurobert

New boundary conditions for the approximate flux-limited diffusion radiative transfer in circumstellar environments Test case study for spherically symmetric envelopes

Keywords:

Context. In order to constrain the models describing circumstellar environments, it is necessary to solve the radiative transfer equation in the presence of absorption and scattering, coupled with the equation for radiative equilibrium. However, solving this problem requires much CPU time, which makes the use of automatic minimisation procedures to characterise these environments challenging. Aims. In this context, the use of approximate methods is of primary interest. One promising candidate method is the flux-limited diffusion (FLD), which recasts the radiative transfer problem into a non-linear diffusion equation. One important aspect for the accuracy of the method lies in the implementation of appropriate boundary conditions (BCs). We present new BCs for the FLD approximation in circumstellar environments that we apply here to spherically symmetric envelopes. Methods. At the inner boundary, the entering flux (coming from the star and from the envelope itself) may be written in the FLD formalism and provides us with an adequate BC. At the free outer boundary, we used the FLD formalism to constrain the ratio of the mean radiation intensity over the emerging flux. In both cases we derived non-linear mixed BCs relating the surface values of the mean specific intensity and its gradient. We implemented these conditions and compared the results with previous benchmarks and the results of a Monte Carlo radiative transfer code. A comparison with results derived from BCs that were previously proposed in other contexts is presented as well.

Results. For all the tested cases, the average relative difference with the benchmark results is below 2% for the temperature profile and below 6% for the corresponding spectral energy distribution or the emerging flux. We point out that the FLD method together with the new outer BC also allows us to derive an approximation for the emerging flux. This feature avoids additional formal solutions for the radiative transfer equation in a set of rays (ray-tracing computations). Conclusions. The FLD approximation together with the proposed new BCs performs well and captures the main physical properties of the radiative equilibrium in spherical circumstellar envelopes.

Introduction

The study of circumstellar environments at different stages of stellar evolution is of crucial importance. These environments reflect the physical processes in action, from the star formation with the presence of accretion discs to late stages in the evolution, in which strong stellar winds shape the circumstellar envelopes. Observations at high angular resolution allow to probe and characterise the circumstellar material by determining densities, temperatures, abundances, velocity fields, etc. The exploitation of instruments such as the Multi AperTure mid-Infrared SpectroScopic Experiment1 (MATISSE), operating in the mid-IR, or the Atacama Large Millimeter Array2 (ALMA) in the sub-millimetric domain offer complementary views of these environments, that give access to regions close to the central star up to the outer regions using a multi-wavelength approach.

Circumstellar matter is generally composed of a mixture of gas and dust particles that absorbs and scatters the incident stellar radiation. The envelope is then heated by the radiation, and a radiative equilibrium can be reached in which the envelope also emits radiation in the infrared domain.

In order to constrain the models describing circumstellar environments, it is necessary to solve the radiative transfer equation under the assumption of radiative equilibrium. Several numerical techniques exist to solve this problem in one, two, and three dimensions. The Monte Carlo method is popular because it can be adapted to any geometry and can handle many physical processes (see Steinacker et al. 2013 for a thorough review). However, solving the radiative transfer problem requires much CPU time, which makes the use of any automatic minimisation procedure to characterise these environments challenging.

In this context, the use of approximate methods is of primary interest. One promising candidate is the flux-limited diffusion (FLD), introduced by Levermore & Pomraning (1981) (L&P hereafter). This description numerically simplifies the problem by recasting the radiative transfer equation into a non-linear diffusion equation for the mean specific intensity of the radiation field (see Sect. 2).

Physically, the boundary conditions (BCs) for the radiative transfer equation are obtained from the known specific intensity incident on the surface of the object. However, the outgoing in-tensity is not known a priori and has to be obtained from the solution of the radiative transfer problem. It is thus not obvious to find a BC for the mean specific intensity at the surface of the object. A consequent theoretical work has already been done in finding satisfying BCs for the FLD method (Pomraning 1986(Pomraning , 1988)). These BCs were derived with the assumption that a boundary layer could be defined, which might not be true in astrophysical applications where the media may not be seen as an infinite half-space by the radiation for some frequencies. Furthermore, as far as we know, they have never been numerically tested in an astrophysical context.

The FLD approximation has already been implemented in several astrophysical applications. Sonnhalter et al. (1995); Yorke & Sonnhalter (2002) used the FLD to solve the frequencydependent radiative transfer to model protostellar discs and massive star formation, respectively. In these studies, the central star was treated as an additional source and the specific mean intensity J ν at the outer edge of the media was set to be equal to the Planck function B ν (T out ), with a prescribed temperature T out at the boundary. Some improvements were made later, in the context of the radiation hydrodynamics problem for massive star formation (see Kuiper et al. 2010;Mignon-Risse et al. 2020). These more sophisticated hybrid codes split the radiation field into two components, the stellar and the dust component, where the FLD method only solves for the latest part. In this treatment, the Dirichlet boundary condition at the outer edge only applies to the dust component. This relies on the assumption that the dust temperature is known at the interface with the interstellar medium. In the problem we consider, the temperature at the outer boundary is not known a priori and must be derived as part of the solution to the radiative transfer problem coupled with the radiative equilibrium equation. In a non-grey problem, we need BCs that can properly handle several regimes of optical thicknesses, for different frequencies.

We present new BCs for the FLD theory in circumstellar environments that we tested and implemented in the case of spherically symmetric envelopes. The condition is derived from the prescription of (i) the incident flux, derived from an extended stellar source and the self-heating of the envelope at the inner boundary and from (ii) the ratio of the mean specific intensity over the radiative flux at the free outer boundary. We show that they both may be written as mixed BCs relating the mean specific intensity and its gradient at the surfaces of the envelope. They consequently lead to a more realistic description of the radiation field (compared to simple Dirichlet or von Neumann boundary conditions) while still remaining sufficiently easy to implement. As a by-product of our investigations, we also derived an approximate expression for the emergent flux at the free outer surface.

The paper is organised as follows: in Sect. 2 we recall the bases of the FLD theory. In Sect. 3 we present the new boundary conditions and in Sect. 4 their numerical implementations. In Sect. 5 we test the accuracy of our results by comparing the temperature profile in the envelope and the spectral energy distribution (SED) of the outgoing flux with the results of two radiative transfer codes, namely (i) DUSTY (Ivezic & Elitzur 1997), which numerically solves the integral equation for the energy density, and (ii) a Monte Carlo (MC) radiative transfer code (Niccolini & Alcolea 2006). Additionally, we compare the derived boundary conditions with the original boundary conditions of Levermore & Pomraning (1981). Finally, in Sect. 6, we conclude and present some perspective for our future work.

The Flux Limited Diffusion Theory

In the following, we present the original work of L&P and introduce the relevant background for the derivation of the BCs in Sect. 3. The position is denoted by r, the direction of propagation by n, and the frequency by the subscript ν. The transport equation for the specific intensity I ν (r, n, t) at the position r in the n direction with isotropic and coherent scattering is

1 c ∂ t I ν + n.∇I ν = -κ ext ν I ν + κ abs ν B ν + κ sca ν J ν .
(1)

J ν = J ν (r, t) is the mean specific intensity, B ν = B ν (T (r, t))
is the Planck function and κ ext ν , κ abs ν and κ sca ν are the extinction, absorption and scattering coefficients, respectively. The zeroth and first moments of the specific intensity, namely J ν and H ν , are defined as

J ν = 1 4π 4π I ν d 2 n , H ν = 1 4π 4π I ν n d 2 n , (2)
where the integration is performed over all directions. These quantities are linked by the zeroth moment of Eq. ( 1),

1 c ∂ t J ν + ∇.H ν = κ abs ν (B ν -J ν ) . (3)
The FLD approximation is a closure of the system of the moment equations by expressing H ν as a function of J ν . This is done by expressing the specific intensity I ν as a function of the mean specific intensity J ν ,

I ν = J ν ψ ν (r, n, t) , H ν = J ν h ν (r, t) , 1 4π 4π ψ ν d 2 n = 1 , (4)
where h ν is the normalised flux and is expressed as

h ν (r, t) = 1 4π 4π n ψ ν (r, n, t) d 2 n . (5)
The ψ ν function is called the normalised intensity and quantifies the anisotropy of the radiation field I ν . In the optically thin and thick limits, this function reduces to a Dirac distribution and a constant, respectively. The FLD approximation consists of assuming that the anisotropy of the radiation field is a conserved quantity, yielding the expression for ψ ν ,

ψ ν (r, n, t) = 1 1 + (h ν -n) .R ν , (6)
with

R ν (r, t) = -∇J ν ω ν κ ext ν J ν , ω ν = κ abs ν B ν + κ sca ν J ν κ ext ν J ν . ( 7
)
The quantity denoted by R ν plays a key role in the description of the local radiation field in the medium. It expresses the ratio of the (effective) mean free path over the characteristic length of the variation of the mean specific intensity. Consequently, the limits R ν 1 and R ν 1 correspond to the optical thin and thick regimes, respectively. The quantity ω ν is called the effective albedo. It is equal to unity in the absence of true absorption (κ abs ν = 0). L&P showed that, using Eq. ( 6) in Eq. ( 5), h ν is proportional to R ν ,

h ν = λ (R ν ) R ν , λ (R ν ) = 1 R ν 1 tanh R ν - 1 R ν , (8)
where R ν is the norm of R ν and λ (R ν ) is the 'flux-limiting' parameter. Finally, using Eq. ( 8) in Eq. ( 4) allows Eq. ( 3) to be rewritten as a non-linear diffusion equation for the mean specific intensity

1 c ∂ t J ν -∇. (D ν ∇J ν ) = κ abs ν (B ν -J ν ) , (9)
with the non-linear diffusion coefficient

D ν = λ (R ν ) ω ν κ ext ν . ( 10
)
We note that in the FLD approach, the radiative net flux H ν is related to the gradient of the mean specific intensity J ν by H ν = -D ν ∇J ν . It shows some similarities with the Fick law, which applies in the stellar interior, but here the diffusion coefficient depends on a non-linear way on the mean specific intensity. In the optically thick regime, the FLD approximation reduces to a linear diffusion equation, whereas in the optically thin regime, it reduces to an advection equation, as expected.

Boundary conditions

In this section, we specify the time-independent BCs that were implemented for the FLD Eq. ( 9). They are defined at the specific boundaries of a circumstellar shell, namely an inner spherical cavity that is illuminated by an enclosed star, and an outer boundary with no incoming radiation. The BCs for the radiative transfer equation Eq. ( 1) would be obtained by specifying the value of the incident specific intensity I ν on the considered surface. In the framework of the FLD approximation, we face two problems: (i) The FLD equation applies to the mean specific intensity, whereas the physical constraint is on the ingoing specific intensity at the boundaries, and (ii) we also explained in Sect. 2 that the FLD approach implies a specific angular dependence of I ν , given by the function ψ ν (see Eq. ( 6)). This specific dependence would not be consistent with any arbitrary value of the incident radiation field at the surface. Consequently, the actual BCs for the radiative transfer equation are in general incompatible with the FLD solution.

This is expected because the FLD approximation is in principle not valid close to the boundaries of objects. Pomraning (1988) derived the BCs for the FLD equation from the decomposition of the radiative transfer problem into an interior problem described by the FLD equation and an additional boundary layer term. The match of the interior and boundary layer solutions yields the BCs for the FLD equation. A relation of the surface values of the mean specific intensity and its gradient is obtained, but with rather involved coefficients depending on integrals of the Chandrasekhar H function and on the surface value of R ν . In the literature, Dirichlet BCs are often used with a prescribed value for the temperature at the outer boundary of the medium. As already pointed out, the surface temperature is not known a priori and must be derived from the solution of the radiative transfer problem.

Here, we propose to impose BCs for the zeroth and the first angular moments of the radiation field, in a form that is compatible with the FLD approach, in order to ensure a smooth match with the interior FLD solution. We obtained mixed Robin-type BCs that relate the surface values of the mean specific intensity and its gradient, but the coefficients are quite simple analytical functions of the surface value of R ν .

Inner spherical cavity

Expression of the incident flux in the FLD formalism

One example of physical significance is to write a condition on the flux F in ν entering a boundary surface at the position r S ,

F in ν (r S ) = J ν (r S ) ŝ.n≤0 ŝ.n ψ ν (r S , n) d 2 n . (11)
Here, ŝ is the outward normal vector to the surface of the envelope ( ŝ = -r at the inner edge). The right-hand side (RHS) of Eq. ( 11) is the expression of an incoming flux in the FLD formalism. Eq. ( 11) is the general form of the inner BC without any assumption on the geometry of the problem. The integral can be performed analytically using Eq. ( 6) for ψ ν if we assume that the vector R ν is normal to the surface (R ν = ± R ν r), which is exact for spherically symmetric problems. The equation can then be rewritten as a non-linear mixed BC,

F in ν (r S ) = π α ν J ν (r S ) + β ν ŝ.∇J ν | r=r S ( 12
)
with

α ν = 2 ln (cosh R ν ) R ν tanh R ν r=r S , β ν = 2D ν | r=r S . (13)
This BC may be regarded as an implicit relation of the surface values of J ν and R ν that yields an FLD solution compatible with the given incident flux. It is interesting to note that although Eq. ( 12) is different in its coefficients and construction from the original BC Eqs. ( 56) and ( 66) in L&P, they are analytically equivalent. In L&P, the coefficients were derived to give an exact transport result for the case of a source-free (B ν = 0) half-space media, with a constant κ ext ν , ω ν and a particular incident intensity distribution Γ(r S , µ) of the form

Γ(r S , µ) = 1 coth(R ν (r S )) -µ = R ν (r S ) ψ ν (r S , µ) . ( 14
)
This specific form is proportional to the angular dependence in the FLD formalism. The correspondence between α ν in Eq. ( 13) and Eq. ( 56) in L&P is

γ = α ν -2λ(R ν ) ŝ.R ν 2 -4λ(R ν ) ŝ.R ν . (15)
Hence, specifying the FLD flux at the boundaries of a spherically symmetric domain will actually give the same BC as is obtained by solving the exact transport result of a source-free (B ν = 0) half-space media, with a constant κ ext ν , ω ν and a particular incident intensity distribution Γ(r S , µ), given by Eq. ( 14). We note that this equivalence only applies for the spherical and planar symmetric systems and no longer holds when we specify Eq. ( 11) for other configurations, where the radiative flux is not orthogonal to the boundary surface.

Incident flux from an extended source and the envelope

We want to write the BC at an inner spherical cavity located at a distance R in from the centre of a star of radius R (r S = R in r) and surface temperature T . For this, we need to specify the flux F in ν . We have two contributions, The first contribution, denoted by B ν (T ) comes from the star and the other from the inner boundary itself and is expressed as

F in ν = ŝ.n≤0 ŝ.n B ν (T ) + I e ν (r, -n) d 2 n . (16)
I e ν (r, -n) = J ν (r) ψ ν (r, -n).
As shown by Fig. 1, the vector r = r S -2(r S .n) n corresponds to the opposite point at the inner boundary, along n. Because of this dependence, this BC is no longer local by nature and cannot be expressed in a closed form, except in spherical symmetry where J ν (r) = J ν (r S ) and ψ ν (r, -n) = ψ ν (r S , -n). To perform the angular integration, we aligned the n z axis with the unitary vector r. For the star, the integration on µ = cos(θ) (θ being the angle between n z and n) spans from µ 0 = 1 -(R /R in ) 2 to 1, and for the inner cavity, it spans from 0 to µ 0 ,

F in ν = 2π            1 µ 0 µ B ν (T ) dµ + J ν (r S ) µ 0 0 µ ψ ν (r S , -µ) dµ            . ( 17
)
The incident flux F in ν is then expressed as,

F in ν = π        R R in 2 B ν (T ) + γ ν J ν (r S )        , (18)
with,

γ ν = 2 µ 0 tanh R ν -ln (1 + µ 0 tanh R ν ) R ν tanh R ν r=r S . ( 19
)
The inner BC for the FLD equation in the case of an inner spherical cavity enclosing a star is

(α ν -γ ν ) J ν (R in ) + β ν ŝ.∇J ν | r=R in = R R in 2 B ν (T ) , (20)
with α ν and β ν , the quantities defined by Eq. ( 13).

We can analytically write the FLD solution with this inner boundary condition in the limit where the envelope is optically thin (R ν 1) and compare it with the known analytical solution for the dilution of the mean specific intensity in free space,

J ν = B ν (T ) 2           1 -1 - R r 2           . ( 21
)
In the optically thin limit, the FLD Eq. ( 9) and the BC Eq. ( 20) reduce to

r 2 J ν = const = R 2 in J ν (R in ) , J ν (R in ) = 1 4 R R in 2 B ν (T ) , (22)
respectively. The solution of the FLD equation, in this limit is then

J ν = 1 4 R r 2 B ν (T ) , (23)
in agreement with Eq. ( 21) when (R /r) 2 1. The relative difference in temperature between Eq. ( 21) and the FLD solution Eq. ( 23) is ≈ 15 % at the star surface, 1 % at r ≈ 2.5 R . For r > ∼ 10 R (values for our test cases presented in Sect. 5), this difference becomes negligible ( < ∼ 0.06 %).

Outer boundary

For the inner cavity, we would like to impose that no incoming radiation enters the external shell of the envelope (F in ν = 0). However, in the FLD formalism, the angular dependence of the radiation is given by the specific form of the ψ ν function (Eq. ( 6)). The incoming radiation vanishes only when R ν becomes infinite, which also results in a sharp-peaked distribution for the emerging specific intensity. This is not physically realistic, and this inconsistency is expected because the FLD method is rigorously not valid close to the surface of the object. Another approach is then required to describe the behaviour of the radiation on the external edge. Inspired by L&P and Pomraning (1986), we seek a BC in the form of a closure relation between the mean specific intensity and the radiation flux at the surface, that is,

J ν (R out ) -ζ ν ŝ.H ν (R out ) = 0 , (24)
with ζ ν , a coefficient we need to determine. At the outer edge and without incoming radiation, ζ ν has to be understood as the ratio of the energy density over the emerging flux. This ratio can be expressed in spherical symmetry as

ζ ν = 1 0 I ν dµ 1 0 µ I ν dµ . ( 25
)
It depends on the anisotropy of the emergent radiation field. In the optically thin limit, the radiation field is along the ŝ direction (spherical symmetry), thus I ν ∝ δ(µ -1) and hence ζ ν -→ 1. In the diffusion regime where the emergent field is isotropic, ζ ν -→ 2. In the original study of L&P, ζ ν was chosen to be equal to 2, which means that it correctly describes the optically thick cases where R ν 1. If the BC is to correctly describe different optical regimes, we need ζ ν to be a function of R ν . In the framework of the FLD approximation, I ν = J ν ψ ν and the specific angular dependence of I ν , given by ψ ν , is used to compute the surface value of ζ ν ,

ζ ν = ζ(R ν ) = 2 + α(R ν ) tanh(R ν ) α(R ν ) + 2λ(R ν )R ν . ( 26
)
In the two limits (R ν 1 and R ν 1), we recover ζ ν -→ 1 and ζ ν -→ 2. We note that the latter limit reduces to the boundary condition Eq. ( 56) in L&P. Now we need to specify the value of R ν at the external boundary of the envelope. The behaviour of the radiation at this interface is only dictated by the interior solution because there is no incoming radiation. Consequently, to ensure a smooth match of the BC and the interior FLD solution, we used a second-order extrapolated value for R ν (R out ) = R e ν (see Eq. 41). Because we prescribe a value for R ν at the external edge, R e ν is also used to compute the non-linear diffusion coefficient

D ν (R out ) = D ν (R e ν ) = λ(R e ν )
/ω ν κ ext ν and the coefficients of the numerical scheme Eq.( 34) on the external edge. The outer BC without incoming radiation that we implemented is then

J ν (R out ) + ζ(R e ν ) D ν (R e ν ) ŝ. ∇J ν | r=R out = 0 , (27)
where we have expressed H ν (R out ) = -D ν (R e ν ) ∇ J ν | r=R out . In Sect. 5.2 we compare this BC in an astrophysical application with respect to the original BCs of L&P.

Approximation for the emergent flux

We used an extrapolation of the non-linear diffusion coefficient D ν (R e ν ) to relate the flux at the external edge to the gradient of the mean specific intensity. As there is no incident flux, this provides us with an approximate expression for the emergent flux, given by

F out ν = 4πH ν (R out ) = -4π D ν (R e ν ) ∇J ν | r=R out . ( 28
)
This approximation is tested in Sect. 5. It reproduces the results of different radiative transfer codes that solve the full transfer equation very well, as shown in Fig. 2 Fig. 3.

Radiative equilibrium and warming of the stellar surface

Because of the geometric extension of the star (see Fig. 1), there is part of the radiation that emerges from the envelope that falls back onto the star,

F fall ν (R in ) = 2πJ ν (R in ) -µ 0 -1 µ ψ ν (R in , µ) dµ = 2πJ ν (R in ) R ν tanh R ν ln 1 + tanh R ν 1 + µ 0 tanh R ν -(1 -µ 0 ) tanh R ν , (29)
with R ν being evaluated at r = R in . We used the same conventions as in Sect. 3.1 to perform the angle integration. This part of the flux is hidden from the rest of the envelope, and the radiative equilibrium inside the cavity leads to a warming of the stellar surface (Niccolini et al. 2003). This effect can be quite dramatic, and reach up to 30 % of the total stellar luminosity that is obscured for optically-thick grey shell, as in the test case presented in Sect. 5.2. This has to be taken into account to properly ensure the radiative equilibrium condition throughout the full space from the stellar surface to the outer boundary of the envelope. To fulfil the radiative equilibrium condition at the stellar surface we write

σT 4 eff = σT 4 + R in R 2 ∞ 0 F fall ν (R in ) dν . ( 30
)
We imposed a fixed value for T eff , and the temperature of the star T was updated accordingly.

Numerical implementation

The FLD Eq. ( 9) is a non-linear diffusion equation that has to be solved numerically for each point of space, time, and frequency.

In the following, we limit ourselves to the 1D time-independent FLD equation,

1 r 2 ∂ r r 2 D ν ∂ r J ν = κ abs ν (J ν -B ν ) . (31)
Here, r denotes the radial variable from the centre of the envelope. An additional constraint to Eq. ( 31) is given by the equation of the radiative equilibrium

∞ 0 κ abs ν B ν dν = ∞ 0 κ abs ν J ν dν . ( 32
)
These equations are solved for a spherically symmetric envelope of inner radius R in and outer radius R out , surrounding a star of radius R .

Numerical scheme

Finite-difference approximation

Following a finite-difference procedure, we discretise and sample in a logarithmic way the frequency domain into n ν points, denoted by the subscript k. Space is discretised into n x cells and denoted by the subscript i. J ν is computed at the cell centres and the vector H ν on the walls. The differential operator is approximated with a second-order finite differences operator. The equation is solved with respect to a new variable x = f (r) on a regular grid of constant step ∆x = (x(R out )x(R in )) / (n x -2) to allow r to be non-uniformly sampled. We make use of one ghost cell for each grid border to ensure the BCs. We obtain the following system of equations:

A k,i+ 1 2 J k,i+1 -A k,i J k,i + A k,i-1 2 J k,i-1 = -b k,i B k (T i ) . (33)
The non-linear coefficients A are given by

A k,i± 1 2 = r 2 i± 1 2 D k,i± 1 2 dx dr i± 1 2 , A k,i = A k,i+ 1 2 + b k,i + A k,i-1 2 , b k,i = ∆x 2 dx dr i -1 r 2 i κ abs k,i . (34)
The non-linear nature of the equation arises from the expression of the coefficients A and RHS in Eq. ( 33). They implicitly depend on J ν , through the diffusion coefficient D ν and the radiative equilibrium Eq. ( 32), respectively. The coefficients A require an estimation of J ν and its gradient (see Eqs. ( 10) and ( 7)) at the cell walls, given by

J k,i+ 1 2 = 1 2 J k,i+1 + J k,i , ∇J ν | k,i+ 1 2 = dx dr i+ 1 2 J k,i+1 -J k,i ∆x r . (35)

Iterative scheme

Several strategies are possible in order to solve the FLD Eq. ( 33) coupled with Eq. ( 32). The simplest approach is to use an iterative method to fully solve Eq. ( 33) and to update the temperature through Eq. ( 32). Iterating between these two processes until convergence yields the solution of the problem. This procedure, commonly called the Λ-iteration in the literature (see Yorke & Sonnhalter 2002), becomes very slow and does not converge for large optical depths. In analogy with the usual accelerated Λ-iteration (ALI) methods, Yorke & Sonnhalter (2002) found an improved convergence behaviour by splitting the solution of equation Eq. ( 31) in Eq. ( 32).

We found a simple method, inspired by the Gauss-Seidel approach, to solve Eqs. ( 33) and ( 32) simultaneously instead of repetitively. If we denote by n the iteration of the method, we have

J n+1 k,i = b k,i B k (T n i ) + A n k,i+ 1 2 J n k,i+1 + A n k,i-1 2 J n+1 k,i-1 A n k,i , (36)
and the temperature is updated after only one Gauss-Seidel spatial sweep at each frequency through Eq. ( 32), which we rewrite as

n ν -1 k=0 W k κ abs k,i B k (T n+1 i ) = n ν -1 k=0 W k κ abs k,i J n+1 k,i . (37)
We replaced the frequency integration by a quadrature formula with the associated weights W k . The left-hand side (LHS) of Eq. ( 37) was pre-computed and stored in a table, for a wide range of temperatures, allowing the RHS to be linearly interpolated in this table (in the logarithm of the integral for better accuracy)

. By doing so, we avoided using a Newton-Raphson procedure to determine the new temperature, which reduces the computational time.

Our procedure consisted of repeatedly updating J ν and T with the help of Eqs. ( 36) and ( 37) until we reached convergence. The coefficients A Eq. ( 34) and the BCs Eq. ( 40) were immediately updated for each frequency k after one Gauss-Seidel spatial sweep. We note that this procedure is different from the usual Λ-iteration presented above because the temperature is updated simultaneously with J ν , within the same iteration n.

Update of the stellar temperature

The radiative equilibrium inside the inner cavity requires updating the stellar temperature (see Sect. 3.4). Following Eq. ( 30), we updated the stellar temperature at the end of each iteration n of the numerical scheme,

T n+1 =         T 4 eff - 1 σ R in R 2 n ν -1 k=0 W k F fall,n+1 k         1 4 . ( 38
)
Here again, we replaced the frequency integration by a quadrature formula with the associated weights W k . F fall,n+1 k is given by Eq. ( 29) and computed with the freshly updated values of J n+1 k, 1 2 and R n+1 k, 1 2 .

Boundary conditions

We used two ghosts cells (one at each boundary of the domain) in order to simplify the implementation of our BCs. In doing so, the inner BC was imposed at the wall between the first (i = 0) and second cell (i = 1), and at the outer BC between the last (i = n x -1) and penultimate cell (i = n x -2). As indicated by Eqs. ( 20) and ( 27), we need to specify J ν and ∇J ν at these interfaces. For this, we used Eq. ( 35) and write

J ν (R in |R out ) ≈ J k,0|n x -2 + J k,1|n x -1 2 , ŝ.∇J ν | R in |R out ≈ dx dr 1 2 |n x -3 2 J k,1|n x -1 -J k,0|n x -2 ∆x . ( 39
)
Table 1. Results from the comparison with DUSTY τ ν 0 (T ) (λF λ /F) p = 0 p = 2 p = 0 p = 2 1 0/0/1 0/0/1 1/1/4 1/1/3 10 0/0/1 1/1/3 2/3/14 2/2/10 100 1/0/1 1/1/3 1/1/8 1/2/11 1 000 2/1/4 1/1/4 6/8/30 3/2/9

Notes. Relative differences (in %) for the temperature profiles (T ) and for the SEDs (λF λ /F) shown in Fig. 2. The results are presented in the form mean( ) / std( ) / max( ) and rounded to the closest percent.

Accordingly, the values of J ν in the ghosts cells were updated immediately after one Gauss-Seidel sweep Eq. ( 36) to ensure the BCs,

J k,0 = 2∆x R R in 2 B k (T ) -∆x (α k -γ k ) -2 dx dr | 1 2 β k J k,1 ∆x (α k -γ k ) + 2 dx dr | 1 2 β k , J k,n x -1 = 2 dx dr | n x -3 2 ζ(R e k )D k (R e k ) -∆x J k,n x -2 2 dx dr | n x -3 2 ζ(R e k )D k (R e k ) + ∆x , (40)
with α k , β k and γ k being defined by Eqs. ( 13) and ( 19). For the extrapolated value R e k in ζ(R e k ) Eq. ( 26) and D k (R e k ), we used a second-order Lagrange extrapolation,

R e k = 3 R k,n x -5 2 -R k,n x -7 2 + R k,n x -9 2 . (41)

Initial conditions

Initial conditions for both J ν and T must be provided in order to solve Eq. ( 33). It is clear that the overall convergence speed strongly depends on the initial setup of the solution, but there is also a trade-off with the stability, that is, the ability of the solution to converge, for a wide variety of cases. As an initial guess, we used the analytic solution of the FLD in the optically thin limit Eq. ( 23), and we write

J 0 k,i = 1 4 R r i 2 B k (T ) , (42)
from which we deduce the corresponding temperature profile T 0 with the help of Eq. (37).

Numerical tests: spherically symmetric envelopes

5.1. Benchmarks from Ivezic et al. (1997) We tested the accuracy of our FLD code with our Robin-type mixed boundary conditions in a general and realistic case, by comparing it with the 1D benchmark problems realised by Ivezic et al. (1997). We recall the conditions of the test and refer to the original paper for further information.

A point source surrounded by a spherically symmetric envelope of matter at radiative equilibrium irradiates as a black body at the temperature T = 2 500 K. This envelope extends from the inner radius R in to the outer radius R out = 1 000 R in . The inner radius is set so that the temperature at the inner radius is always T in = T (R in ) = 800 K. The density profile n(r) is assumed to be

10 0 10 1 10 2 10 3 r /R in 10 -1 10 0 T /T in 10 0 10 1 10 2 10 3 r /R in 10 -1 10 0 10 1 10 2 10 3 10 4 λ [µm] 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 λF λ /F 10 -1 10 0 10 1 10 2 10 3 10 4 λ [µm]
Fig. 2. Non-grey case : Normalised temperature profiles (upper panels) and SEDs (lower panels) for four different opacities τ ν 0 = 1, 10, 100, and 1 000 (blue, orange, green and red, respectively) and two density power laws: p = 0 (left panels) and p = 2 (right panels). The solid lines represent the FLD curves, and the black dots indicate the benchmark profiles from Ivezic et al. (1997).

a power law of the form n(r) = n 0 (R in /r) p . The radial optical depth τ ν of the envelope is linked to the density profile by

τ ν = R out R in κ ext ν dr = R out R in C ext ν n(r) dr , (43)
where C ext ν is the extinction cross-section coefficient. It is defined by

C abs ν = C abs ν 0 , C sca ν = C sca ν 0 if ν ≥ ν 0 , C abs ν = C abs ν 0 ν ν 0 , C sca ν = C sca ν 0 ν ν 0 4 if ν ≤ ν 0 , C ext ν = C sca ν + C abs ν , (44)
with

C abs ν 0 = (1 -η) C ext ν 0 , C sca ν 0 = η C ext ν 0 ,
ν 0 the frequency corresponding to λ = 1 µm and η the albedo, set to 1/2 for these tests. The benchmark problems are thus completely defined by two parameters: (i) the exponent in the density power law p = 0, 2, and (ii) the radial optical depth of the envelope at ν 0 , τ ν 0 = 1, 10, 100, and 1 000. This created eight different cases to test the accuracy of our code. The coefficient n 0 in the density profile is derived with the help of τ ν 0 and p,

n 0 = (p -1) τ ν 0 C ext ν 0 R out R out R in p        R out R in p-1 -1        -1 . (45)
The normalised temperature profile T/T in and the normalised SED λF λ /F (F = ∞ 0 F λ dλ) of the envelope are shown in Fig. 2 for each case.

The Ivezic benchmarks were produced with version 2 of DUSTY3 (Ivezic & Elitzur 1997). Because our code is of different nature, the spatial and frequency grids are different. We then compared the results by linearly interpolating our profiles (in loglog scale) on the DUSTY grids. We used 128 points for space and frequency, with a logarithmic sampling. The corresponding relative differences are displayed in Table 1. We also point out that we restricted the comparison, for the normalised SEDs, to the frequency domain where λF λ /F ≥ 10 -6 because of non-physical results of the DUSTY code below this threshold, for the smallest wavelengths.

The FLD results and the benchmarks agree well. The average relative differences in the temperature profiles mean( (T )) is of the order of 1 %, with a maximum value of approximately 4 %, achieved by the most optically thick envelopes (τ ν 0 = 1 000). The average of the relatives differences in the normalised SEDs mean( (λF λ /F)) always stays below 3 %, with the exception of the optically thick envelope with constant density profile, where this difference reaches 6 %. 5.2. Grey spherical shell with the Monte Carlo code from Niccolini & Alcolea (2006) We compared the FLD code with a 3D MC radiative transfer code (Niccolini & Alcolea 2006). We wished to test the new BCs Levermore & Pomraning (1981) (dashed lines) for the test case presented in Sect. 5.2. The colour code is the same as in Fig. 3. The emerging fluxes are displayed in scale to highlight the differences between the BCs. The lower panels show the relative differences profiles with respect to the MC code from Niccolini & Alcolea (2006) in temperature (left) and in emerging flux (right).

in the less realistic but more extreme case of a spatially small spherically symmetric grey envelope. We expect the boundary effects to play a major role for this type of problems. The inner radius was set to R in = 10 R and the outer radius R out = 20 R . We assumed a constant density profile (p = 0) in the envelope. Our test cases consisted of determining the normalised tempera-ture profiles T/T in and the emerging fluxes F λ for several cases, ranging from optically thin (τ = 0.01) up to the optically thick envelopes (τ = 100). The corresponding profiles are shown in Fig. 3.

Because the codes are different, we interpolated our results linearly (in loglog scale) on the MC grids. The relative differ- (T ) (F λ ) 0.01 0/0/0 0/3/17 1 1/1/5 0/3/17 10 1/1/7 1/8/33 100 0/0/3 1/2/13 Notes. Relative differences (in %) for the temperature profile (T ) and for emerging flux (F λ ) shown in Fig. 3. The differences are presented in the form mean( ) / std( ) / max( ) and rounded to the closest percent. mean( ) is given by Eq. ( 46). ences between the two codes, are displayed in Table 2. As an additional feature, the MC code also provides an estimation of the errors on the temperature σ(T ) and on the emerging flux σ(F λ ), computed from the MC noise (Niccolini & Alcolea 2006). We used this information to compute a more relevant mean value for (T ) and (F λ ),

mean( (T )) = N x i=0 W i (T i ) N x i=0 W i , mean( (F λ )) = N λ k=0 W k (F k ) N λ k=0 W k (46)
where N x (N λ ) is the number of spatial (wavelength) points of the MC grid, (T i ) ( (F k )) is the relative error (in %) on the temperature (emerging flux) between our results and the MC results, and W i (W k ) is the inverse square of the MC relative errors, defined as

W i = σ(T i ) T i -2 , W k = σ(F k ) F k -2 (47)
The two results agree well. The average of the relative differences of the temperature profile mean( (T )), remains of the order of 1 % for all the cases we tested. The largest differences are reached for the intermediate cases (τ = 1, 10) and are located on the external edge of the envelope. This is expected because the FLD approximation is known to perform well in the optically thin and thick regimes, but it is less well suited to describe these intermediate cases. Nevertheless, the temperature profile is still quite well reproduced and the emerging flux is not affected by the small errors on the temperature close to the outer edge. We point out that the BCs derived in this paper allow us to successfully reproduce the correct behaviour of the temperature profile for the optically thick envelope where it shows a quite steep decrease at the outer surface, as shown by Fig. 3. The emerging fluxes F λ agree within 1 % on average, except for the optically thick case (τ = 100), where mean( (F λ )) reaches about 5 %. We conclude this section with a comparison of the outer BC from this study Eq. ( 27) and the original BC Eq. ( 56) from L&P, given by,

J ν (R out ) + 2D ν ŝ.∇J ν | r=R out = 0 (48)
We previously mentioned in Sect. 3.1 that the inner BC Eq. ( 20) is analytically identical to that of L&P in spherical symmetry, so we restrict the comparison to the outer edge of the envelope.

It is important to notice that the BC of L&P was not originally intended to describe this class of problems, but the comparison still remains instructive for studying the importance of the BCs for the accuracy of the solution. In Eq. ( 48), the factor 2 (Eq. 58 in L&P) was originally used to give the correct ratio of the energy density over the emerging radiative flux in plane-parallel geometry for an optically thin slab illuminated by an isotropic incident radiation field. However, for the special case of a spherical envelope surrounding a black-body star, this ratio becomes

1 µ 0 B ν (T )dµ 1 µ 0 µB ν (T )dµ = 2 (1 -µ 0 ) 1 -µ 2 0 = 2 1 + µ 0 (49) with µ 0 = 1 -R R out 2
, the cosine of the stellar angular size at R out . As µ 0 → 0 or equivalently R /R out → 1, this ratio increases to 2, as expressed by L&P, because we recover the case of a plane-parallel geometry with an isotropic incident radiation field. Far from the star (µ 0 → 1 or R /R out → 0), the ratio tends to 1, associated with an incoming sharp-peaked radiation. Hence, the L&P BC that set the ratio to 2 will strongly deviate from the analytic limit 1, in the optically thin regime. In Fig. 4 we display the relative differences in the temperature profiles (T ) and in the emerging fluxes (F λ ) for the same test case as presented at the beginning of this section. We note that this test case is not realistic, however, it allows to compare different optical regimes, in contrast to a more realistic problem in which the radiation is free in the external regions most of the time, such as for the test case presented in Sect. 5.1. We recall that the BC from L&P is a limiting case of the BC derived in this study, in the case R 1, where the emerging radiation is almost isotropic (see Sect. 3.2), hence it is not surprising that the results converge to the same profile, for a optically thick grey envelope (τ = 100). On the other-hand, in the optically thin case (τ = 0.01) where we would expect the ζ coefficient Eq. ( 26) to be close to unity, the BC from L&P performs poorly as expected. For intermediates regimes (τ = 1, 10), the entanglement of the error of the BC and FLD method itself makes any comparison very hard. We note that although the temperature profile is closer to the benchmark result for τ = 10 with the BC of L&P, this is not the case for the associated emerging flux. We also note that the BC, although defined locally, can have a global effect on the whole solution, as shown by the temperature profile of the test case τ = 1. To conclude, we also point out that we tried to implement the L&P BC in the non-grey cases (Sect. 5.1) but we were unable to reach a satisfying convergence of the computations.

Conclusion

The FLD approximation together with the new BCs, yields promising results in correctly describing the radiation transport inside spherically symmetric circumstellar envelopes. These conditions, derived in Sect. 3 from physically consistent constraints on the behaviour of the radiation field at the inner and outer surfaces, allow us to compute with a good accuracy the temperature profile and the SED for a wide range of configurations, from very small to very large optical thicknesses. As shown in Sect. 5, it reproduces the correct temperature profile within ≤ 2%, and the SED or emerging flux at less than ≤ 6% on average, with respect to the solution of the full radiative transfer equation under radiative equilibrium.

The numerical solution of the 1D non-linear diffusion equation Eq. ( 31) coupled with the radiative equilibrium Eq. ( 32) was performed with a Gauss-Seidel method-based iterative scheme, in which the temperature is updated at each iteration step. Furthermore, we point out that the FLD approximation implemented with the proposed outer-boundary condition provides a simple approximation for the flux emitted by the envelope. This allows computing the SED or emerging flux without using any ray-tracing module, which is a significant gain in computational time.

The next step, which will be our main concern for our future work, is the generalisation of these boundary conditions to non-spherically symmetric media, in particular for circumstellar discs, where the BCs will take a more complex form. In this regard, the extension to 2D will require numerical optimisations. Several directions of improvements are already being studied, such as the use of multi-frequency adaptive spatial grids, or acceleration procedures for efficiently solving the FLD equation under the radiative equilibrium condition. The last point is also crucial for solving extreme optically thick envelopes.

Fig. 1 .

 1 Fig. 1. Geometry of an inner spherical cavity (solid black line) illuminated by a central star (grey disc).

Fig. 3 .

 3 Fig. 3. Grey case : Normalised temperature profiles (left panel) and emerging fluxes (right panel) for four different opacities τ = 0.01, 1, 10, and 100 (blue, orange, green and red, respectively) with a constant density profile (p = 0). The solid lines represent the FLD curves, and the black dots with the error bars σ indicate the MC profiles from Niccolini & Alcolea (2006).

Fig. 4 .

 4 Fig. 4. Comparison of the outer boundary conditions presented here (solid lines) and fromLevermore & Pomraning (1981) (dashed lines) for the test case presented in Sect. 5.2. The colour code is the same as in Fig.3. The emerging fluxes are displayed in scale to highlight the differences between the BCs. The lower panels show the relative differences profiles with respect to the MC code fromNiccolini & Alcolea (2006) in temperature (left) and in emerging flux (right).

  τ

Table 2 .

 2 Results from the comparison with the MC code

http://www.eso.org/public/teles-instr/paranal-observatory/vlt/vltinstr/matisse

2 http://www.almaobservatory.org

available at http://faculty.washington.edu/ivezic/dusty_web/

Acknowledgements. This study has been supported by the Lagrange laboratory of Astrophysics and funded under a 3-year PhD grant from Ecole Doctorale Sciences Fondamentales et Appliquées (EDSFA) of the Université Côte-d'Azur (UCA). The benchmark profiles were generated with the radiative transfer code DUSTY, available at http://faculty.washington.edu/ivezic/dusty_web. Part of the computations were carried out with the help of OPAL-Meso computing facilities. The authors are grateful to the OPAL infrastructure from Observatoire de la Côte d'Azur (CRIMSON) for providing resources and support.


Asef Ghabeli
email: asef.ghabeli@savba.sk

Enric Pardo
email: enric.pardo@savba.sk

Mark Ainslie
email: mark.ainslie@eng.cam.ac.uk

Loïc Quéval
email: loic.queval@geeps.centralesupelec.fr

Modeling HTS dynamo-type flux pumps: open-circuit mode and charge of an HTS coil

Keywords: Dynamo-type flux pump, open-circuit DC voltage, full-circuit mode, REBCO high temperatuere superconductor

Dynamo-type flux pumps can be employed to energize high-temperature superconducting (HTS) magnets to maintain the steady current mode or in rotor windings of electrical HTS motors and generators to inject the current without the need for current leads. To understand the complicated operating working mechanism of the flux pump it is necessary to study both in opencircuit and connected to an HTS coil. In this work, we employ two models to analyze these configurations. The first one is the Minimum Electromagnetic Entropy Production (MEMEP). The second one is a segregated H-formulation finite-element model. For open-circuit mode, we investigate the effect of the airgap on DC voltage generation. Later, we also study the energizing process of an HTS coil modeled by a series resistance and an inductance. These models will also be compared to experiments.

I. INTRODUCTION

HTS dynamo-type flux pumps belong to the category of the traveling magnetic wave type, which uses a moving magnet to create an inhomogeneous traveling magnetic flux across an HTS tape in order to generate voltage. These kinds of flux pumps are simple and easy to maintain because they do not utilize any switches. For this reason, they have become very popular during recent years to be exploited in HTS magnets, motors and generators [START_REF] Walsh | Characterization of current stability in an hts nmr system energized by an hts flux pump[END_REF]- [START_REF] Bumby | Development of a brushless hts exciter for a 10 kw hts synchronous generator[END_REF]. Dynamo-type flux pumps can be employed to energize HTS magnets to maintain the steady current mode or in rotor windings of electrical HTS motors or generators to inject the current without using any brushes in a contactless way.

Among the different types of flux pumps, the dynamo-type is one of the most complicated ones in terms of understanding its mechanism. Modeling has showed to be an effective way to explain the flux pump mechanism. For this reason, several articles have been published in literature to explain the dynamo-type flux pump mechanism [START_REF] Campbell | A finite element calculation of flux pumping[END_REF]- [START_REF] Ghabeli | Modeling of airgap influence on DC voltage generation in a dynamo-type flux pump[END_REF].

This work is funded from the Slovak VEGA agency, under contract No. 2/0097/18.

To investigate the flux pump mechanism, studying both open-circuit and full-circuit mode is essential. The open-circuit mode is simpler compared to full-circuit mode considering the elimination of the effect of a net current flowing in the HTS circuit. In this work, we utilize mainly the Minimum Electromagnetic Entropy Production (MEMEP) [START_REF] Pardo | Electromagnetic modelling of superconductors with a smooth current-voltage relation: variational principle and coils from a few turns to large magnets[END_REF], [START_REF] Pardo | 3d computation of non-linear eddy currents: Variational method and superconducting cubic bulk[END_REF] to model a dynamo-type flux pump with realistic superconductor properties (E -J power-law relation with J c (B, θ) from critical-current measurements of single tapes) [START_REF] Ghabeli | Modeling of airgap influence on DC voltage generation in a dynamo-type flux pump[END_REF]. This variational method solves the local current density as the state variable, only inside the superconducting region, with minimization of a functional, which makes the model fast and efficient. In addition, the results are compared with another one implemented in COMSOL based on a segregated Hformulation finite-element model [START_REF] Quéval | Superconducting magnetic bearings simulation using an h-formulation finite element model[END_REF]. This model solves the magnetic field strength H only inside the superconductor and a thin air layer around it, which allows us to reduce the number of degrees of freedom and the computing time.

This work investigates the effect of airgap on DC voltage generation in a dynamo-type flux pump in open-circuit mode.

In addition, we study the charging procedure by the dynamotype flux pump of an inductance and resistance in series, representing a typical REBCO coil with joint resistance.

II. AIRGAP DEPENDENCE OF OPEN-CIRCUIT VOLTAGE

Fig. 1 shows an overview of the modeling configuration. The flux pump rotates in the xy-plane and the magnetic field B exist in the same plane. The current density J flows in the z direction. Therefore, the vector potential A and electric field E are parallel to the z axis.

A useful and practical quantity measured experimentally is the output voltage difference between the tape at 77 K (superconducting mode) and at 300 K (normal mode), ∆V = V 77K -V 300K . This quantity ∆V contains the DC component of the open-circuit voltage, ∆V oc . The amplitude of ∆V oc decreases with increasing airgap, which is the result of the reduction in magnetic field density on the tape surface. Fig. 2 shows the calculated results of reduction of ∆V oc for the modeled flux pump from 2.4 mm up to 50 mm. By calculating the ∆V oc signals in Fig. 2 for the airgaps ranging from 2.4 to 50 mm, the stationary average DC voltage per cycle is

III. CHARGING PROCESS OF HTS COIL

Studying the full-circuit mode with transport current is essential to understand the flux pump mechanism that energizes superconducting windings. Fig. 4 shows the circuit diagram of the flux pump in open-circuit and connected to an HTS coil, respectively. In Fig. 4(left) the net transport current is zero, and hence the generated voltage by the flux pump V F P is equal to V oc as calculated in the previous section. However, in the full circuit of Fig. 4(right), the output voltage for the flux pump is related to the transport current by

V F P (I(t), t) = RI(t) + L coil dI(t) dt (1)
where R denotes the joint resistance and L coil is the load inductance. Solving the non-linear differential equation above requires numerical methods. Using similar arguments as in [START_REF] Li | Coupling loss at the end connections of rebco stacks: 2d modelling and measurement[END_REF], the full circuit mode in MEMEP can be solved by minimizing the functional in [START_REF] Pardo | Electromagnetic modelling of superconductors with a smooth current-voltage relation: variational principle and coils from a few turns to large magnets[END_REF] for the current density at each time step with the additional term RI 2 /2+L coil (∆I) 2 /2, where ∆I and ∆t are the time and current increments between two time steps, and ∆I is found during minimization. For the segragated H-formulation model, ( 1) is solved together with the H-formulation flux pump model using a bidirectional circuit coupling scheme. Modeling results will also be compared to experiments.

Fig. 1 .

 1 Fig. 1. This 2D cross-sectional model configuration is sufficient to describe the main features of a dynamo-type flux pump.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Calculated results showing the trend in the reduction of ∆Voc with increasing airgap up to 50 mm for the case with Jc(B, θ) dependence and rotating frequency of 12.3 Hz. (∆Voc as calulated in [7])..

Fig. 4 .

 4 Fig. 4. Circuit diagram of a dynamo-type flux pump in open-circuit mode (left) and in full-circuit mode (right). Voc is the open-circuit voltage and R denotes the joint resistance. L coil denotes the load inductance.


domain: spi phys

text:
Quentin Herbette
email: quentin.herbette@onera.fr

Stéphane Saillant

Michel Menelle

Bruno Urbani

Gilbert Auffray

Nicolas Bourey

Muriel Darces

Marc Hélier

HF Radar antenna near field assessment using a UAV

Keywords: Unmanned Aerial Vehicle, Near Field, HF Surface Wave Radar

This paper deals with a system aboard a UAV used to map the electric near field, in amplitude and phase, radiated by an antenna operating in the HF band (3-30 MHz). The objective is to characterize the radiation pattern of antennas that are commonly integrated in HF surface wave radar. An innovative and original approach is required because usual characterization methods do not suit well for that frequency band. The proposed device is described in the paper and associated constraints are highlighted. A trial of that system has been performed on a monopole antenna, and comparisons between measurements and simulations are in good agreement.

I. INTRODUCTION

ONERA is currently deploying biconical antennas for its sky wave and surface wave radars. A photography of a part of the transmitting array of a HF radar operated by ONERA, composed of numerous biconical antennas, is shown below in Fig. 1. The large size of HF antennas, and difficulties to model the environment in which they operate, complicate the applications of conventional antenna characterization techniques, as anechoic chamber or Open Area Test Site. Even if HF antenna far field characterization by means of a system aboard a helicopter has already been achieved, that kind of tests are very expensive and difficult to perform [START_REF] Stehle | RELEDOP: a full-scale antenna pattern measurement system for high-powered HF transmitting antennas[END_REF]. To address this problem, a new method and a new airlifted system are proposed. It will allow to obtain the radiation pattern of an HF antenna by means of a near-field to far-field transformation, thanks to a previously developed software [START_REF] Djoma | Prediction of Sky and Surface Wave Radiation of a Wideband HF Antenna[END_REF]. One of the features of that software is its ability to take the surface wave radiation into account, which is required to characterize HF antennas, used for surface wave radar, in their own environment. In the future, as the proposed device is flexible, it will also be applied to new innovative transmission systems that are still in development. In section II, the principle of measurement and the main constraints that apply to the proposed device are presented. In section III, the solution to the synchronization problem between the transmitter and receiver is explained. In section IV, the final system placed aboard the UAV is described. And finally, comparisons between the measured radiation pattern of a monopole antenna and a simulated one are commented in section V.

II. PRINCIPLE OF MEASUREMENT AND ASSOCIATED

CONSTRAINTS

To correctly characterize the near field of an antenna, it is required to map the three components 𝐸 𝜌 , 𝐸 𝜑 and 𝐸 𝑧 of the near electric field, in amplitude and phase . This information must be known for different geographical points located on the surface of a virtual volume centered on the antenna as it can be seen in Fig. 2. Fig. 2. Illustration of the measurement surface centered on the transmitting system composed by three biconical antennas Similar measurement has already been achieved in the laboratory for higher frequency (1 GHz) by applying the electromagnetic equivalence principle, using an electric near field probe named EFS-105 designed by Enprobe [START_REF] Bourey | In Situ Antenna Far Field Estimation Based on Equivalent Sources[END_REF].

In order to achieve the near-field to far-field transformation previously described [START_REF] Djoma | Prediction of Sky and Surface Wave Radiation of a Wideband HF Antenna[END_REF], it is imperative to quantify the electric field phase. It is therefore necessary, by some kind of synchronization, to estimate the phase difference between the antenna under test feed signal and the received signal. However, at 10 MHz, a precision on phase shift of 10° corresponds to a precision of 3 ns on time shift. This issue is the major constraint for the development of the near-field measurement system.

Another critical issue is the precision needed for the near field probe positioning, in order to minimize the error induced on the calculated far-field. By assessing the robustness of the developed software at 1 GHz, we were able to determine that the needed precision was at least ± λ 30 (±1 m at 10 MHz). This matter was quickly resolved by integrating a GNSS RTK (Global Navigation Satellite System Real Time Kinematic) positioning system into the UAV. It provides an accuracy of less than 10 cm on the drone positioning. Finally, the use of a UAV imposes two other restrictions. First, the weight of the components needs to be optimized to reduce the impact on the drone autonomy. Last, it is necessary to move the electric field probe further apart from the drone in order to mitigate coupling effects between the carrier and the near field probe.

III. SYNCRHONIZATION

Performing a synchronization of the digitization clocks of the transmitted and received signals is a major challenge which is yet to be resolved [START_REF] Orraca | Compact high frequency receiver system for Arecibo observatory radiation pattern measurement[END_REF]. Indeed, a nanosecond precision is difficult to achieve, even with a miniature atomic clock or with a GPSDO (Global Positioning System Disciplined Oscillator). That issue is related to the inevitable clock drift that affects the operating frequency of those tools as well as their sensitivity to vibrations, which are not negligible while using an UAV.

Therefore, it has been necessary to explore an alternative allowing to acquire the electric field phase without an equipment such as the one described above.

The chosen solution is to bring back, at ground level, the signal acquired by the near field probe, which is aboard the drone, by means of an optical fiber. This method has been concurrently developed by other research teams [START_REF] Fritzel | Introduction into an UAV-based near-field system for in-situ and large-scale antenna measurements[END_REF] and it allows us, by integrating an electro-optical converter on the UAV, to gather near electric field values, without digitizing the signal aboard the drone. Moreover, the optical fiber is not sensitive to radio frequency radiation and therefore resilient to any disturbances on the electric field information. At last, lightweight optical fiber has a low impact on the weight of the payload. This design makes possible to mimic the manipulation performed in laboratory and simplifies data processing. The next step has been to design the payload.

IV. PAYLOAD CONCEPTION

The UAV used is a DJI Matrice 600 pro. This kind of industrial drone can lift as much as 6 kg for 15 minutes fly duration. The payload is mounted under the drone and it is necessary to comply with the integration restrictions in weight and size. To this end, an electro-optical converter has been modified in order to reduce its size and weight, so that it can be integrated inside the payload. For optimizing the payload weight, it has been decided to choose aluminum to build the platform hosting the different parts of the system. The payload, Fig. 3, is made up of the electro-optical converter (①), a 12 V battery (②) and the near field probe case housing (③). To mitigate coupling effects between the transmitting antennas and the drone, it is necessary to move the near field probe aside from the UAV. The electromagnetic simulation software CST (Computer Simulation Technology) has been used to estimate the needed distance between the sensor and the drone. Unfortunately, the drone manufacturer does not communicate on the properties of the UAV materials. Consequently the simulation has been performed assuming the worst-case scenario, which is that of a perfect electrical conductor (PEC) structure. In order to have a negligible error, it is required to locate the sensor at least at 70 cm from the UAV. For practical reasons due to landing constraints, the probe has to be placed on one of the sides of the drone. The required distance is obtained with a fiberglass tube as shown in Fig. 4. The last part designed is the clamping device used to position the sensor in the three directions corresponding to the three electric field components. In order to simplify the measuring system and avoid adding extra weight, it has been decided not to use a remotely controllable motor that can change the orientation of the probe. The operation is therefore manually performed on ground. The drone positioning accuracy is about ten centimeters, so it can be considered that the UAV will stand at quasi-unchanged measurements points for the three flights. The clamping device, Fig. 5, has been manufactured in polycarbonate.

A. Measurement setup

In order to check the proper functioning of the abovedescribed system, it has been chosen to acquire the radiation pattern of a quarter-wave monopole antenna. This well-known kind of antennas enables to perform simple comparisons between measurements and simulations. A 6 m monopole antenna has been erected on a 20×25 m² ground plane. This ground plane was mounted above a concrete ground with unknown electrical parameters. The antenna impedance was matched to 50 ohms at a frequency of 13.5 MHz thanks to a RLC matching circuit. To avoid damages to the near field probe, which is in its saturation region for an electric field higher than 30 V/m, the antenna input power has been limited to 6 dBm. Feed and received signals have been both digitized with a NI PXIe-5171 oscilloscope that allows a 14-bit quantization and a 250 MS/s sample rate. Both signals have been digitized using the same clock, so the relative phase difference between the two signals is known. Measured antenna and associated ground plane can be seen in Fig. 6.

Considering the radiation pattern of the monopole antenna as omnidirectional, it was not needed to measure the near electric field on a cylinder centered on the antenna. Indeed, by only measuring it on a radial line and duplicating the obtained data, it has been possible to generate a closed measurement cylinder around the antenna. GPS flight data have been recovered after each flight and have been used to determinate the position the probe at any time.

B. Comparisons between measured and simulated near electric field

The three electric field components have been measured. Since the near electric field 𝐸 𝜑 associated to a monopole antenna is close zero, results are given only for 𝐸 𝜌 and 𝐸 𝑧 electric field components. As shown in Fig. 7, two trajectories have been followed to measure 𝐸 𝜌 and 𝐸 𝑧 . After neglecting slight movements in y plane, comparisons with simulations have been performed.

CST Microwave and FEKO have been used to simulate near electric field radiated by the monopole antenna exhibited in figure 6. The antenna reflection coefficient 𝑆 11 was measured in order to know the exact radiated power. In consequence, comparisons between measurements and simulations have been made with absolute amplitude data.

Electrical delay, due to the coaxial cable and optical-fiber lengths, is not precisely known, so the field phase comparisons are achieved by shifting the measured phase with an arbitrary constant electrical delay. So, the differences between measured and simulated phases are interesting to compare. In addition, absolute field phase values are not needed to perform the near field to far field transformation with the previously developed software. Horizontal axis represents the total distance flown by the UAV. Zero corresponds to the position of the antenna. The first 25 m coincide with the displacement along the x-axis by the drone at a constant altitude. The last 9 m corresponds to the drone descent phase and landing. Comparisons between measured and simulated 𝐸 𝑧 data (amplitude and phase values) show a good agreement. However, there are some differences between the simulated data from CST and FEKO: that might be explained by the different electrical properties of the concrete used in simulation. A considerable difference is shown on a part of the 𝐸 𝜌 phase. It has not yet been elucidated, but it is of little importance. Indeed, the associated amplitude values are small and will have a low impact on the near-field to far-field transformation. Several successive measurements were made and the obtained results are similar. Near electric field measurements were satisfactory and reproducible, and so it has been possible to achieve a near field to far field transformation.

C. Near field to far field transformation and comparisons

The obtained far field radiation pattern has been compared to the simulated one computed using FEKO. Far field has been evaluated at 300 m from the monopole antenna. As it can be seen on figure 10, there is a good agreement between the two radiation patterns. In both cases, the maximum radiated power is at 𝜃 = 50°. Differences between the two curves can be explained by the lack of knowledge about the properties of the concrete covering the ground on which the antenna is located. Those results confirm the relevance of carrying out a characterization of HF antennas with the proposed system. Successful measurements on a monopole operating at 50 MHz were also performed and we look forward to test our measuring system on VHF array to determine the upper operating frequency of the system. In this paper, a drone-embedded measuring system has been proposed, designed and tested. It is used to acquire the near electric field radiated by an antenna operating in the HF band. In order to perform the near-field to far-field transformation with a previously developed software, the amplitude and phase of the three components of the near electric field are required. The innovative solution here is to overcome the difficulty of synchronizing the clocks of transmission and reception systems, mounted aboard the UAV, by recovering the field values via an optical fiber. The constraints that apply to the design of the payload have been reviewed and the devices fulfilling them have been assembled and installed in the payload. Last, a successful test of the entire system has been carried out on an actual HF monopole antenna.

Fig. 1 .

 1 Fig. 1. Biconical antennas of a HF radar transmitting system operated by ONERA

Fig. 3 .

 3 Fig. 3. Payload placed under the UAV

Fig. 4 .

 4 Fig. 4. DJI Matrice 600 pro equipped with the measurement system

Fig. 5 .

 5 Fig. 5. Clamping device used to position the probe V. COMPARISONS BETWEEN MEASUREMENTS AND SIMULATIONS

Fig. 6 .

 6 Fig. 6. Measured monopole antenna and associated ground plane

Fig. 7 .

 7 Fig. 7. Actual drone flight paths respectively used for measurements and simulations of 𝐸 𝜌 and 𝐸 𝑧

Fig. 9 .

 9 Fig. 9. Comparisons between measured and simulated near electric field 𝐸 𝜌 , in absolute amplitudes and relative phases

Fig. 10 .

 10 Fig. 10. Comparisons between measured and simulated far field evaluated at 300 m from the antenna VI. CONCLUSION

ACKNOWLEDGMENT DGA (French Ministry of Defense) is thanked for its contractual support to the RAPID AC3M project for which the airborne system has been designed and developed.