Agnes Jocher

Heinz Pitsch

Thomas Gomez

Jérôme Bonnety

Guillaume Legros
email: guillaume.legros@upmc.fr

Combustion instability mitigation by magnetic fields

Keywords: numbers: 47.15.Fe, 47.70.Pq, 47.65.Cb

come

I. INTRODUCTION

Spontaneous low frequency oscillations of atmospheric non-premixed coflow flames were found to be induced by buoyancy-driven Kelvin-Helmholtz (KH) instabilities [1]. These instabilities can be triggered by a shear layer in a flow and might appear in the atmosphere, oceans, and stellar internal flows [2]. The prediction and the subsequent control of stability limits is of significant interest not only in flames and nature, but also in a wide range of industrial applications [3], for example the oxyfuel combustion technique [4]. This technique reduces anthropogenic CO 2 release into the atmosphere due to high O 2 concentrations and recycled combustion products -mainly water vapor and CO 2 -in the unburnt oxidizing stream. Doing so, the peak soot content generated along the combustion process is also decreased. A major drawback of the technology is the enhanced tendency to combustion instability, especially caused by CO 2 addition. This leads to limited combustion reliability and in the worst case to mechanical failure and damage of the whole system. Furthermore, in unstable combustion regimes, the likelihood of unburnt fuel and soot expulsion increases due to local flame quenching and suppressed oxidation processes. The harmful effects of soot particles on human health and environment have been discussed extensively [5][6][7] and have led to the tightening of regulations by governments on emission sources like vehicles and power plants as well as on the ambient concentration itself [START_REF] Agency | Environmental Protection Agency -Control of Emissions from Nonroad Spark-Ignition Engines and Equipment; Final Rule[END_REF][START_REF] Commission | Communication from the Commission to the European Parliament, the Council, the European Economic and Social Committee and the Committee of the Regions -A Clean Air Programme for Europe[END_REF][START_REF] Agency | Environmental Protection Agency -Control of Hazardous Air Pollutants from Mobile Sources; Final Rule[END_REF]. To help overcoming these drawbacks, we suggest to combine technologies. Although special attention has been paid to the combustion enhancing role of electric fields [START_REF] Hilliard | [END_REF]12] the stability of flames experiencing magnetic perturbations has not been thoroughly examined. It was shown in the literature that the onset of the KH instabilities in non-premixed flames can be suppressed by increased soot production and the subsequent flame cooling due to the associated enhanced radiative emission by soot particles [13]. An upward gradient of the square of the magnetic flux density ∇(B 2 ) applied to a steady laminar non-premixed sooting coflow flame was shown to enhance soot formation in the flame, while the formed soot was still fully oxidized, i.e. no soot was released through the flame tip [14]. Furthermore, experiments [15] evidenced that buoyancy induced convection in non-premixed flames can be partly substituted, enhanced, or reduced by magnetic fields. In several patents, significant reductions in NO x and CO emissions have been documented when permanent magnets were located around the fuel injector of an internal combustion engine [START_REF] Twardzik | Apparatus for subjecting hydrocarbonbased fuels to intensified magnetic fields for increasing fuel burning efficiency[END_REF][START_REF] Menzell | Emission control devices[END_REF][START_REF] De John | Assembly and process for improving combustion emissions of a combustion apparatus[END_REF]. However, no experimental evidence for a nonpremixed flame stabilization and soot emission modification has been shown on an academic configuration up to now.

Here, we document for the first time the experimental observation of stabilizing a spontaneously oscillating non-premixed sooting flame with a magnetic field. The context of oxyfuel combustion technique is addressed. We analyze the observation with the help of numerical simulations, applying the aforementioned magnetic effects on both mass transfer and soot production to reveal an original phenomenology that potentially allows control-ling the stability domain of reacting flows, which is then evidenced by a local stability analysis. The new methodology that could be imagined from the results presented here is especially relevant to oxygen enhanced and carbon dioxide diluted combustion strategies, therefore to low-impact combustion devices.

II. EXPERIMENTAL OBSERVATIONS AND PROCEDURE

The present study was initiated after we had observed an unexpected response to an external magnetic perturbation of a non-premixed flame that exhibits natural low frequency (12.6 Hz) flickering. Indeed, this spontaneous instability vanished as the flame was submitted to a magnetic field generated by an electromagnet. The movie provided as supplementary material [19] displays the time-history of the raw data documenting the initial experiment. The evolution of the non-premixed flame with time is captured by a Phantom v711 highspeed camera equipped with a widescreen CMOS sensor. To record this direct visualization of the experiment, the camera just substituted the one shown in Fig. 1(a). A sequence of 48 s at 188 frames per second of 12-bit monochrome frames is captured on an 800 x 304 pixels 2 matrix. The exposure time was kept constant at 5.3 ms and the camera was focused on the plane containing the flame's vertical axis of symmetry, using a SIGMA 105 mm F2.8 Macro lens. With this optical arrangement, each pixel in the CMOS array focused light from a volume corresponding to 0.13 mm in height, 0.13 mm in width and 1 mm in depth. The voltage V delivered by a photodiode evidences the magnitude of the flames radiative signature along the experiment. Starting with the spontaneously flickering flame when no magnetic field is applied, the current I flowing through the coils of the electromagnet is gradually increased up to 60 A, generating a maximum magnetic field magnitude of 1 T and a maximum ∇(B 2 ) of 18.2 T 2 /m. The voltage evolution delivered by the photodiode evidences that the flame flickering is fully suppressed for the maximum current flowing through the coils of the electromagnet. In the following, experimental, numerical, and theoretical tools are deployed to assess this original stabilization process.

Figure 1(a) displays the setup that enables this experiment together with the associated investigations. As extensively outlined by Jocher et al. [14], the burner used (1) is inserted between both coils (2) of the electromagnet, with the burner's exit plane located 130 mm below the coils' horizontal axis of symmetry. Thus, within the volume occupied by the non-premixed flame investigated here, a constant upward ∇(B 2 ) [14] can be generated with a decently uniform maximum magnitude of 18.2 T 2 /m. In the following, the case without magnetic field is called MagF 0 and the one with |∇(B 2 )|=18.2 T 2 /m is referred to as MagF 1. None of the rig components is magnetic to avoid mechanical in-terferences. By virtue of its design and size, the whole arrangement cannot be considered a practical combustion device as the experimental flame height does not exceed 5 cm, while the outer diameter of the coils is 80 cm (see Fig. 1(b)). Practical combustion devices exhibiting the same effect could be designed smaller than the current academic experiment. However, for appropriate modeling and analysis this academic arrangement was chosen to establish well controlled and spatially constant conditions in terms of ∇(B 2 ).

For the present study, the axisymmetric laminar, flickering, non-premixed ethylene flame shown in the movie burns in 55% oxygen mixed with 45% carbon dioxide in volume, on a Santoro type coflow burner [START_REF] Santoro | [END_REF]. Ethylene is chosen as a fuel because its sooting propensity is documented in a large body of literature and an extensively studied detailed mechanism including polycyclic aromatic hydrocarbon (PAH) chemistry is available [21,22]. The coflow burner consists of concentric brass tubes with effective inner diameters of D f =11 mm and D ox =102 mm. The flow rates of the axial ethylene stream and the coflowing oxidizer are adjusted by two Bronkhorst EL-FLOW mass flow controllers to 0.36 l/min and 74 l/min, respectively. The corresponding Reynolds numbers for these conditions are 87 and 1438. Even though the flow is laminar, the non-premixed flame exhibits spontaneous flickering [23]. It is worth noting that the stabilization process documented hereafter could not be observed for a low sooting, flickering, nonpremixed methane flame at similar fluid dynamic conditions.

To track the soot volume fraction distribution, the Laser Extinction Method (LEM) [23,24] is applied, providing two-dimensional soot volume fraction fields with both fine temporal and spatial resolutions. The arrangement of the optical diagnostics to conduct LEM specifically through a flame located in the electromagnet is shown in Fig. 1(a) and (b). The system consists of a 100 mW continuous wave laser (3) operating at 645 nm (-5/+7 nm) as the monochromatic light source, a neutral density filter (4) to adjust the intensity of the laser, a set of beam expanding optics (6) including a rotating diffusive disc (7) to generate a uniform beam intensity, a set of collection optics [START_REF] Agency | Environmental Protection Agency -Control of Emissions from Nonroad Spark-Ignition Engines and Equipment; Final Rule[END_REF], and a camera [START_REF] Commission | Communication from the Commission to the European Parliament, the Council, the European Economic and Social Committee and the Committee of the Regions -A Clean Air Programme for Europe[END_REF]. A digital pulse generator controls the occurrence and the duration of the CMOS exposure, together with the shutter (4) opening. A frame grabber records on a computer the frames captured by the camera, such as the one shown in Fig. 1(a). The Photon Focus MV1 12-bit progressive scan monochrome camera ( 9) is mounted with a conventional lens equipped with a narrow band filter centered at 645 nm (± 2 nm), as well as with a band width, at one half of the transmissivity maximum of 20 nm. With this optical arrangement, the matrix of 1312 x 1082 pixels 2 provides a spatial resolution of 137 µm for the LEM projected data, over the 60 mm diameter area of investigation. For the current study, the frames were recorded at a frame rate of 94 Hz and an exposure time of 5.4 ms. The flame is considered an emitting, absorbing, but non-scattering medium. For the laminar coflow non-premixed ethylene flames studied, the flame radiative spectrum in the visible is governed by the continuum radiation from soot. This is particularly true in the upper part of the visible spectrum. In addition, absorption by soot particles produced in these non-premixed flames is shown to be at least one order of magnitude higher than scattering, especially at large wavelengths in the visible range [25]. In such a configuration, the Radiative Transfer Equation that models the transfer of the radiative intensity can be integrated along an optical pathway, such as the one followed by the collimated laser beam inside the flame in Fig. 1 allows the attenuation, i.e. the difference between both quantities, to be only connected to the spectral absorption coefficient field κ λ (r, z) crossed along the optical pathway leading to the pixel considered. In practice, at every height z i imaged on a line of pixels (see Fig. 1(c)), deconvoluting the attenuation measurements integrated over the line-of-sight leads to a system of linear equations that is solved for (κ λ ) ij at the locations r j along the line. As every set of equations is shown to be ill-conditioned, a Tikhonov regularization is used to stabilize the deconvolution process [23]. The whole field κ λ (r, z) can then be retrieved. The soot volume fraction field f v (r, z) can also be inferred following the Mie theory and assuming that soot particles are in the Rayleigh limit as

f v (r, z) = [λ κ λ (r, z)] / [6 π E(m)] , (1)
where E(m) is a function of the complex refractive index m of soot. Following the methodology of Kashif et al. [23], E(m) is here adjusted to 0.43 to reproduce the peak soot volume fraction measured by Santoro et al. [START_REF] Santoro | [END_REF] at HAB=50 mm for ethylene burning in air. The level of the mean relative uncertainty within a region where the soot volume fraction exceeds 10% of its peak value is then found to be slightly lower than 5%. However, this does not incorporate the uncertainty associated with the refractive index function of soot particles reflected in the ongoing debate about this quantity.

III. NUMERICAL SIMULATION OF THE OBSERVED STABILIZATION PROCESS

In the following paragraph, the two-dimensional, axisymmetric, numerical setup, and the applied methods to reproduce the experimentally observed stabilization process are outlined. The parallelized, finite difference code CIAO solves the Navier-Stokes equations in the low Mach number limit, utilizing spatial and temporal staggering, together with Crank-Nicolson type time advancement [26,27]. The Poisson equation for pressure is solved by the multi-grid HYPRE solver and the scalar equations are solved via the bounded quadratic upstream interpolation for convective kinematics (BQUICK) scheme, based on the QUICK scheme of Leonard [28]. The temperature and species equations are advanced by introducing a symmetric operator split of Strang [29] and the chemistry operator uses a time-implicit backward difference method similar to that implemented in CVODE [30]. The flow is treated as a multi-component mixture, where diffusion is approximated by the Curtiss-Hirschfelder approach [31] together with a correction velocity to account for mass conservation. The second-order Soret diffusion caused by temperature gradients is taken into account for all species, while heating due to viscous dissipation, as well as the second-order Dufour process are neglected. The chemical mechanism [32], containing 47 species and 290 reactions, was shown to predict soot precursors, including PAH chemistry. Gas phase and soot radiation is implemented with a discrete ordinates method (DOM) [33] including CO 2 , CO, and H 2 O. The Hybrid Method of Moments (HMOM) [START_REF] Mueller | Proc. Combust. Inst[END_REF] predicts soot quantities by taking the volume and surface of the soot particles into account. To account for the magnetic field influence, the body force term F i acting on a chemical species i and its resulting drift are added to the momentum, species, and temperature equations [START_REF] Yamada | [END_REF]. The force can be expressed as

F i = 1 2 µ 0 ρ Y i χ i ∇(B 2 ), (2)
where µ 0 = 4π • 10 -7 is the magnetic permeability of vacuum, Y i the mass fraction of species i, and χ i its magnetic susceptibility per unit mass. The latter quantity is given by Curie's law as

χ i = N A g 2 L µ 2 B S i (S i + 1) µ 0 3 k T m i , (3)
where N A =6.022•10 23 1/mol, g L =2 and µ B =9.274•10 -24 J/T are the Avogadro number, the Lande's g-factor and the Bohr magneton, respectively. k=1.38•10 -23 J/K is the Boltzmann constant, T the temperature, and m i the molar mass of species i. The total angular momentum of the species i electron spin is defined as S i . The Lorentz force on ionic species is assumed negligible, due to the absence of ionizable substances like alkali elements [36]. The numerical setup corresponds to the experimental configuration. The burner nozzle extends 14 mm into the computational domain and the coflow duct's exit plane is located 11 mm into the domain. In the following the origin is set to the intersection of the vertical burner axis of symmetry and the fuel nozzle's exit plane. By choosing a numerical domain of 300 mm in radial and 120 mm in axial direction, it is ensured that the results are not affected by the numerical domain boundaries. At the bottom boundary, inflow conditions are used for the fuel and oxidizer flows, surrounded by stagnating air. Symmetry conditions are applied on the centerline, free-slip on the free-stream side and zero gradient conditions on the top boundary. The mesh is non-uniform, and cylindrical with 192 (z) x 240 (r) control volumes. The minimum resolution is 0.02 mm in both axial and radial directions close to the nozzle outlet and in regions of high temperature and soot volume fraction gradients. The inflow temperature and pressure are set to 298 K and 1 atm, respectively. The oscillating non-premixed flame simulation is initialized from the steady solution and evaluated after 10 cycles, corresponding to 0.64 s in real time. The number of cycles needed to ensure a fully periodic state was identified by a start-up transient analysis [37]. The here employed numerical methods and soot formation models have been validated in an earlier paper [27] by comparing a steady ethylene flame burning in a coflow stream of air to color-ratio pyrometry and laser induced incandescence measurements from two different laboratories. The soot formation and oxidation zones have been well predicted as well as the location and the magnitude of the peak soot volume fraction. The impact of flame oscillations on the soot predictions was investigated in the same paper [27] by pulsing the described steady coflow flame at two frequencies. Again, the trends were predicted really well. As observed in the experiments, the computed flames were oscillating at the prescribed pulsing frequencies and the instantaneous peak soot volume fraction was shown to shift from the centerline to the wing of the flame within one oscillation cycle. Also, the expulsion of the soot pocket and the subsequent flame collapse occurred in temporal agreement with the measurements. Here, the oscillation frequency is not induced by an external flow pulsing but by buoyancy and therefore it depends on other flame parameters, such as the local flame temperature. Still the computed and measured oscillation frequencies are very close at 15.6 Hz and 12.6 Hz, respectively. Furthermore, as in [27], the changes of the peak soot volume fraction in the wings of the oscillating flame can be reproduced numerically. of the soot volume fraction field are over predicted in the simulation. However, comparing the present simulations with state-of-the-art simulations [38] of soot production in non-premixed flames, the agreement between experiments and numerical simulations can be considered good. Finally, and most importantly, the original experimental observation we are focusing on in this paper, a magnetic flame instability mitigation could be reproduced computationally. In both, the experiments and the simulation, the non-premixed flame is stabilized with a ∇(B 2 ) magnitude of 18.2 T 2 /m. Based on the simulations' capability of reproducing the experimentally observed natural flame oscillation as well as the discovered flame stabilization with applied magnetic gradient we are confident that the simulations reproduce the experimental results within the required scope to use the numerical simulation to function as the base result for the stability analysis performed later in this paper.

IV. MECHANISMS DRIVING THE NON-PREMIXED FLAME STABILIZATION

After showing that the experimentally observed instability suppression could be reproduced numerically, the potential mechanisms that drive the phenomenon are now assessed. Figure 3 shows the normalized radial profiles for the axial velocity (a), temperature (b), and ethylene (c) and oxygen (d) mass fractions for three non-dimensional heights above the burner z/D f . The profiles are obtained by the corresponding dimensional variables from the two-dimensional numerical simulation. The variables are time averaged over one cycle (64 ms) with a discrete dataset sampled every millisecond. Then, the radial distance, the axial velocity, the temperature, and ethylene and oxygen mass fractions are non-dimensionalized with respect to the fuel nozzle inner radius r f =5.5 mm, the fuel inflow mean velocity U f,in =0.063 m/s, the inlet temperature T in = 298 K, the unity ethylene mass fraction at the inflow, and the coflow mixture's oxygen mass fraction of 0.47, respectively. Three non-dimensional heights above the burner z/D f =0.4, 0.6, and 1.0 are displayed for both magnetic field cases. For both conditions, the peak velocity increases with increasing z/D f due to the buoyant acceleration of the hot combustion gases. The observed changes in the non-dimensional mean flow profiles are consistent with the findings by Jocher et al. [39]. Among the species present in the combustion process, the major contribution to the magnetic force is attributed to the paramagnetic oxygen molecules at ambient conditions in the coflow, due to their high mass fraction at comparably low temperatures [START_REF] Yamada | [END_REF]. The resulting magnetic force per unit volume f mag that initiates and sustains the so-called thermo-magneto convection can then be expressed as [39]

f mag = ρ Y O2 χ O2 -ρ (∞) Y (∞) O2 χ (∞) O2 ∇ B 2 2 µ 0 . (4)
The superscript (∞) indicates here the conditions away from the flame. As shown by Jocher et al. [39], the term between brackets in Eq. ( 4) is always negative for non-premixed flames. As a result, a significant thermomagneto convection can be sustained and is here opposed to buoyancy due the upward direction of ∇ B 2 . The global residence time inside the flame is then increased when increasing the magnitude of the upward ∇(B 2 ) leading to enhanced soot production. Consequently, the thermal expansion is also enhanced when moving from MagF 0 to MagF 1, leading to the outward shift of the velocity, temperature, and mass fraction profiles. The reduced peak velocity combined with the radially outward shift at MagF 1 leads to a reduced peak shear layer as well as shear layer broadening (not shown here), therefore, to a possible suppression of the KH instability. Jocher et al. [14] argued that due to the partial compensation of buoyancy by thermo-magneto convection, the particle residence time in the flame is increased, leading to a higher peak of the radially integrated soot volume fraction. Katta et al. [13] found that the level of soot radiation that is enhanced due to such an increase can induce a significant local flame cooling. Therefore, the buoyant acceleration of hot combustion gases is weakened and the flickering flame stabilized. Combining both findings points to the conclusion that due to the enhanced soot formation and subsequent radiation, the flame temperature decreases and therefore reduces the acceleration of hot combustion gases due to buoyancy. Consequently, the shear layer between the hot combustion gases and the cold coflow is reduced, together with the flame's sensitivity to the KH instability. The latter finding could explain why the KH instability was suppressed by the magnetic gradient in a sooting non-premixed ethylene flame, while no such stabilization could be observed for a low sooting non-premixed methane flame at similar fluid dynamic conditions.

V. LOCAL STABILITY ANALYSIS

Finally, the growth rates of small perturbations in the non-premixed flame at MagF 0 and MagF 1 are investigated by a spatial, local, inviscid stability analysis in cylindrical coordinates, assuming a low Mach number and a parallel and swirl-free mean flow in axial direction. The main goal is to show that this basic theory can qualitatively predict the influence of the magnetic gradient on the spatial growth rates of small perturbations in the flame. When discussing the perturbation's growth rate -α i obtained by a local spatial stability analysis, it is important to keep in mind that it describes the instability behavior of the local mean flow, but not of the entire flow field. Still, it provides a framework to identify the external excitation frequency β that could be necessary to control the flow [40]. An inviscid analysis is performed and justified here based on the findings of Mahalingam et al. [41] that -α i at a finite Reynolds number is usually lower than in the inviscid case. The governing equations for continuity and momentum [42] are complemented by equations for temperature and fuel and oxidizer mass fractions. For the chemical source term, a one-step reaction model is implemented as described by Mahalingam et al. [41]. The influence of the magnetic body force implementation into the governing equations used for the stability analysis was found to be negligible compared with the changes in the computed mean flow profiles. For the same reason, the gravitational force is not considered in the governing equations of the stability analysis [1]. Therefore, magnetic field and buoyancy influences actually impact the stability analysis results only through the mean flow obtained from the two-dimensional simulation, where both effects are considered. The normalized mean flow profiles for the axial velocity U m , the temperature T m , and the ethylene Y C2H4 and oxygen Y O2 mass fractions shown in Fig. 3 are used as input profiles for the stability analysis. Additional input parameters are the molecular weights of ethylene and the oxidizer. The stoichiometric coefficients applied here are 1 and 3, respectively. Similar to Mahalingam et al. [41] the non-dimensional value of the activation energy is kept equal to 3, the heat release parameter to 1000, and the Damköhler number to 8, for all cases shown here. The Lewis and Prandtl numbers are assumed unity and the Reynolds number infinity [41]. In addition, the Soret and Dufour effects, diffusion due to pressure gradients, and the effects of radiation are not considered in the perturbation equations. Due to the low Mach assumption, the constant thermodynamic pressure appears in the thermal equation of state, hence, ρ m is directly proportional to 1/T m . In the following, only axisymmetric perturbations (m=0) are considered. The perturbations are assumed three-dimensional, while the mean flow is one-dimensional in streamwise direction.

The resulting growth rates -α i of this generalized eigenvalue problem of β range are shown in Fig. 4 for the mean flow profiles in Fig. 3. For every height z/D f studied, the peak -α i is higher for MagF 0 as compared to that for MagF 1. This suggests that the flow is more prone to develop an instability without the magnetic gradient exposure. Furthermore, the neutral stability point at which -α i approaches zero is always found at lower β for MagF 1. Thus, the magnetic gradient exposure reduces the unstable frequency range. Although, according to the experiment, the growth rate for MagF 1 should be negative over the whole range of frequencies, it is worth reminding that the conducted analysis addresses local stability and that a viscous stability analysis will show lower amplification rates [41]. For a more quantitative agreement, extensive efforts should be devoted to the development of a global stability analysis. Still, the basic theoretical tool that we contributed here to extend unveils a stabilizing tendency supporting the experimental observations.

VI. CONCLUSIONS

We conclude that the exposure to an upward gradient of the square of the magnetic flux density can reduce a non-premixed flame's sensitivity to KH instabilities by increasing the soot content and soot expulsion to the en-vironment due to flame quenching in such an oscillating non-premixed flame can then be inhibited. Furthermore, we show that the flow modification by a magnetic gradient severely changes the non-premixed flame base structure and consequently soot production in the flame. Reduced growth rates for perturbations in a non-premixed flame under a magnetic gradient exposure could be identified by a local stability analysis. We anticipate our investigation to be a starting point for more sophisticated combustion control techniques. In future investigations the steady non-premixed sooting flame will be submitted to a downward gradient of the square of the magnetic flux density to assess a possible soot formation reduction induced by the resulting magnetic body force and the gravitational force pointing into the same direction. Another potential implementation could be the systematic shaping of magnetic fields [43,44] to achieve a well defined modification of combustion processes.

FIG. 1 .

 1 FIG. 1. (a): Schematic of the arrangement allowing the experiment: (1) coflow burner; (2) coils of the electromagnet; (3) continuous wave laser; (4) shutter; (5) neutral density filter; (6) set of beam expanding optics; (7) rotating diffusive disc; (8) set of collection optics; (9) camera. (b): Close-up view of the laser beam crossing the flame in between both coils. (c): Typical frame captured by the camera as the shutter is open. The burner tip can be seen at the bottom and the soot layer produced in the flame right above.

  (a) and (b). When the laser is off (shutter closed), the energy E (of f ) λ accumulates on a pixel of the camera during a time ∆t due to the steady impinging flux emitted by the flame at wavelength λ. When the laser is on (shutter open), the energy E (on) λ accumulated on the same pixel is E (of f ) λ complemented by the energy deposited by the incident non-coherent light ray. Measuring consecutively E (on) λ and E (of f ) λ

Figure 2

 2 displays in the left-most four images the measured (a) and computed (b) soot volume fraction fields for one oscillation cycle without magnetic field. The local soot volume fraction is normalized by the peak soot volume fraction f v,max of 17 ppm and 63 ppm in the experiments and the numerical simulations, respectively.The peak soot volume fraction and hence the axial extension

FIG. 2 .

 2 FIG. 2. Normalized soot volume fraction fields fv/fv,max. Experiments (a) and simulations (b). The left-most four images show oscillating fields without magnetic field, the isolated image on the right shows the steady flame stabilized by a magnetic field. The phase angle is computed as φ = 360 • f ∆t, where f is the respective frequency (experiment: 12.6 Hz; simulation: 15.6 Hz) and ∆t the physical time after the occurrence of the maximum flame height.

FIG. 3 .

 3 FIG. 3. Normalized computed radial mean flow profiles. (a) Axial velocity Um, (b) temperature Tm, and (c) ethylene YC 2 H 4 ,m and (d) oxygen YO 2 ,m mass fractions at z/D f =0.4 (circle), 0.6 (cross), and 1.0 (triangle) without (solid line) and with (dashed line) exposure to the steady magnetic field gradient.

FIG. 4 .

 4 FIG.4. Perturbation growth rate -αi as a function of the perturbation frequency β. Without (solid line) and with magnetic field influence (dashed line) for z/D f =0.4 (circle), 0.6 (cross), and 1.0 (triangle).

ACKNOWLEDGMENTS

This work was performed using HPC resources from GENCI-CINES (Grant 2015 -c20132b7144). It was financially supported by the Université francoallemande/Deutsch-Französische Hochschule. The authors are grateful to M. Kashif, H. Dutilleul, and J. M. Citerne for their technical support, and to Lutz Lesshafft for providing the base code for the stability analysis and his valuable comments.


Anna E Dudek
email: aedudek@agh.edu.pl

Lukasz Lenart
email: lenartl@uek.krakow.pl

Keywords: almost periodically correlated, autocovariance function, Fourier coefficients, periodicity, testing MSC 2010: 62G09, 62G10, 62M10

In this paper a subsampling approach for nonstationary time series with a non-zero mean function is proposed. It is applied for periodically and almost periodically processes. Two statistical tests are constructed. An example with real data is presented.

Introduction and problem formulation

In this paper we discuss a subsampling approach for nonstationary time series, which have a non-zero mean function. First we introduce some notation. Let {X t , t ∈ Z} be a nonstationary time series and µ(t) = E(X t ) be its mean function. We assume that µ(t) ≡ 0. Moreover, we consider the following decomposition of the mean function µ(t) = µ 1 (t) + µ 2 (t).

Functions µ 1 (t) and µ 2 (t) can be interpreted differently depending on the application. In economics they can represent e.g., seasonal or business fluc-tuations or a long-term trend, while in telecommunications these can be two periodic functions e.g., the first one may be a pilot tone and the latter one the mean function of some periodic signal. A pilot tone is a signal transmitted for control, synchronization or reference purposes. In the following we assume that µ 1 (t) is non-zero and we know how to estimate it. On the other hand, µ 2 (t) may be a zero function and moreover we may not have complete knowledge of its form and properties. Such a model describes a real data situation in which only partial information about the mean function is available e.g., we know the nature of some structure contained in the data like seasonality, but other components require further investigation. We denote the observed sample by X n = (X 1 , X 2 , . . . , X n ). Let μ1,n (t) be the estimator of µ 1 (t) based on X n . We consider the modified sample Xn = ( X1,n , X2,n , . . . , Xn,n ), where Xt,n = X t -μ1,n (t).

In the sequel, for the sake of simplicity we skip the subscript n and we denote elements of the triangular array Xn by Xi , i = 1, . . . , n. From now we will be using only the sample Xn i.e., we consider partially or entirely demeaned data. This corresponds to the standard procedures of data analysis in which the data is assumed to be demeaned before any further analysis (e.g. secondorder) is performed.

In the next section we present how the subsampling idea for nonstationary time series (see chapter 4 in [START_REF] Politis | Subsampling[END_REF]) can be adapted to our setting. In Section 3 we construct two tests designed for time series with periodic and almost periodic structure. The first one is devoted to the firstorder frequency detection, while the second one can be used for testing the second-order periodicity. Finally, Section 4 contains a real data application of our results.

Subsampling

In the sequel we use notation introduced by [START_REF] Politis | Subsampling[END_REF]. Let θ = θ(P) ∈ R be the parameter of interest, where P is the joint probability law governing {X t , t ∈ Z}. Parameter θ can simply relate to µ 2 (t), but one may consider other characteristics of X t . Moreover, let θn = θn ( Xn ) be the estimator of θ and θn,b,t = θn,b,t ( Xt , Xt+1 , . . . , Xt+b-1 ) be the estimator of θ based on subsample Xt , Xt+1 , . . . , Xt+b-1 . In addition, let θ n,b,t = θn,b,t (X t , X t+1 , . . . , X t+b-1 ), where

X i = X i -µ 1 (i), i ∈ Z.
By J b,t (P ) and J n (P ) we denote the sampling distributions of τ b ( θ n,b,t -θ)

and τ n ( θn -θ), where τ (•) is an appropriate normalizing sequence. The corresponding cumulative distribution functions are defined as

J b,t (x, P ) = P τ b ( θ n,b,t -θ) ≤ x (1)
and

J n (x, P ) = P τ n ( θn -θ) ≤ x . (2)
Finally, following [START_REF] Politis | Subsampling[END_REF] we assume that X t is α-mixing.

Conditions below are sufficient conditions for the subsampling consistency.

Assumption 2.1 There exists a limiting law J(P ) such that (i) J n (P ) converges weakly to J(P ) as n -→ ∞.

(ii) For any continuity point x of J(P ) and for any sequences n, b with n, b -→ ∞ and b/n -→ 0, we have

1 n-b+1 n-b+1 t=1
J b,t (x, P ) -→ J(x, P ). OR For any continuity point x of J(•, P ) and for any index sequence {t b } we have that J b,t b (x, P ) -→ J(x, P ) as b -→ ∞.

(iii) max 1≤t≤n-b+1 τ b θn,b,t -θ n,b,t p -→ 0.
One may notice that comparing to [START_REF] Politis | Subsampling[END_REF] (see Assumption 4.2.1. in [START_REF] Politis | Subsampling[END_REF]) we need an additional assumption (iii). It is caused by the fact that in contrast to [START_REF] Politis | Subsampling[END_REF] we are using partially demeaned data. Finally, condition (i) is identical to the one proposed by Politis et al., while (ii) is adjusted to our problem.

The approximation of J n (x, P ) is defined as

Ln,b (x) = 1 n -b + 1 n-b+1 t=1 1 τ b θn,b,t -θn ≤ x ,
where 1{A} is an indicator function of the event A. Below we state the asymptotic validity of subsampling for general statistics.

Theorem 2.1 Assume that {X t , t ∈ Z} is α-mixing time series. Under Assumption 2.1 and taking

τ b /τ n -→ 0, b = b(n) -→ ∞ such that b/n -→ 0 as n -→ ∞, we have (i) if x is a continuity point of J(•, P ), then Ln,b (x) p -→ J(x, P ), (ii) if J(•, P ) is continuous, then sup x∈R | Ln,b (x) -J(x, P )| p -→ 0, (iii) for α ∈ (0, 1), let c(1 -α) = inf{x : J(x, P ) ≥ 1 -α} and cn,b (1 -α) = inf{x : Ln,b (x) ≥ 1-α}. If J(•, P
) is continuous at point c(1-α), then subsampling confidence intervals are asymptotically consistent i.e.,

P τ n ( θn -θ) ≤ cn,b (1 -α) -→ 1 -α. (3)
Remark: In our methodology we applied subsampling to partially demeaned data i.e., to Xn . However, an alternative approach is possible. One may consider the usual subsampling and use subsamples (X t -μ 1,n,b,t (t), X t+1 -μ1,n,b,t (t + 1), . . . , X t+b-1 -μ1,n,b,t (t + b -1)), where μ1,n,b,t (•) is the estimator of the function µ 1 (•) based on (X t , X t+1 , . . . , X t+b-1 ). In such cases the subsampling consistency can be obtained under the sufficient condition presented in [START_REF] Politis | Subsampling[END_REF] (see Assumption 4.2.1). However, one should be aware that the estimator μ1,n,b,t (•) may have higher variance or bias than μ1,n (•), especially when b is small.

Applications

In this section we present possible applications of our subsampling approach. In the first subsection we introduce the class of time series with periodic and almost periodic structure. Additionally, we discuss the assumptions that we use in the consecutive subsections in which we provide two tests.

PC/APC time series

In our consideration we focus on periodically correlated (PC) and almost periodically correlated (APC) time series, which are used for modeling signals in many different fields like economics, telecommunication, mechanics, vibroacoustics and hydrology (see e.g. [START_REF] Hurd | Periodically Correlated Random Sequences: Spectral Theory and Practice[END_REF], [Napolitano, 2012], [Napolitano, 2016]). Time series X t is called PC with period T if it has periodic mean and autocovariance functions, i.e.

µ (t) = µ (t + T ) and B(t, τ ) = Cov (X t , X t+τ ) = Cov (X t+T , X t+τ +T )
for each t, τ ∈ Z. For more details we refer the reader to [START_REF] Hurd | Periodically Correlated Random Sequences: Spectral Theory and Practice[END_REF]. Sometimes it may happen that period length is not known or a considered signal is a composition of two components with incommensurable periods.

To model such data APC time series are used (see e.g. [Napolitano, 2012]). They are a generalization of PC sequences. A time series that has finite second order moments is called APC if its mean and autocovariance functions are almost periodic in time. Finally, a function f (t) : Z -→ R of an integer variable is called almost periodic if for any > 0 there exists an integer L > 0 such that among any L consecutive integers there is an integer p with the property sup t∈Z |f (t + p ) -f (t)| < (see [Corduneanu, 1989] for more details). Analysis of PC and APC processes is often performed in the frequency domain. In such cases the following Fourier expansions of the mean and the autocovariance functions are used

µ(t) = ψ∈Ψ m(ψ)e iψt , B(t, τ ) = λ∈Λτ a(λ, τ )e iλt (4)
for any shift τ ∈ Z, where

m(ψ) = lim n→∞ 1 n n j=1 µ(j)e -iψj , a(λ, τ ) = lim n→∞ 1 n n j=1 B(j, τ )e -iλj and Ψ = {ψ ∈ [0, 2π) : m(ψ) = 0} and Λ τ = {λ ∈ [0, 2π) : a(λ, τ ) = 0}
are the sets of the frequencies. The expansions (4) hold, if the sets Ψ and Λ τ , τ ∈ Z are finite. In the case of PC time series this condition always holds and additionally, we know that Ψ, Λ τ ⊆ {2kπ/T, k = 0, 1, . . . , T -1} for any τ ∈ Z.

Below we present the conditions that are used in the sequel. Let {X t , t ∈ Z} be an α-mixing time series. Depending on the considered case we will assume that X t is PC or APC.

A1 There exists a constant δ > 0 such that

sup t∈Z E |X t | 2+δ < ∞ and ∞ k=0 α(k) δ/(2+δ) < ∞; A2 {X t , t ∈ Z} is uniformly bounded i.e., there exists a constant M > 0 such that P (|X t | ≤ M ) = 1 for all t ∈ Z. Moreover, there exists a constant ζ > 0 such that ∞ k=0 (k + 1) 2 α(k) ζ/(8+ζ) < ∞; A3 For any τ 1 , τ 2 , τ 3 ∈ Z the function (t, τ 1 , τ 2 , τ 3 ) -→ E (X t X t+τ 1 X t+τ 2 X t+τ 3 ) is periodic in t i.e., E (X t X t+τ 1 X t+τ 2 X t+τ 3 ) = E (X t+T X t+τ 1 +T X t+τ 2 +T X t+τ 3 +T ) .

Testing the significance of a frequency in the mean function

In this section we assume that {X t , t ∈ Z} is an APC time series. Let us recall that we consider µ(t) = µ 1 (t) + µ 2 (t) and we know the estimator of µ 1 (t). Moreover, functions µ 1 (t) and µ 2 (t) have the following Fourier expansions µ

1 (t) = ψ∈Ψ 1 m(ψ)e iψt and µ 2 (t) = ψ∈Ψ 2 m(ψ)e iψt , where Ψ = Ψ 1 ∪ Ψ 2 , the set Ψ 1 is known and Ψ 1 ∩ Ψ 2 = ∅.
Then the natural estimator of µ 1 (t) has the form (see e.g. [START_REF] Hurd | Periodically Correlated Random Sequences: Spectral Theory and Practice[END_REF])

μ1,n (t) = ψ∈Ψ 1   1 n n j=1 X j e -iψj   e iψt . ( 5
)
Our aim is to test if for some frequency, for which we know that ψ 0 ∈ Ψ 1 , we have that ψ 0 ∈ Ψ 2 . This problem can be equivalently expressed in terms of the corresponding Fourier coefficients i.e., we consider H 0 : |m(ψ 0 )| = 0 vs H 1 : |m(ψ 0 )| > 0. In this case θ = |m(ψ 0 )| and the test statistics T n based on Xn is of the form

T n = √ n | mn (ψ 0 )| = 1 √ n n j=1 (X j -μ1,n (j)) e -iψ 0 j .
Large values of T n suggest that H 1 is true. The quantiles of T n can be approximated using the subsampling method. For that purpose below we state the subsampling validity for |m(ψ 0 )|.

Theorem 3.1 Assume that condition A1 is fulfilled. Moreover, let the set Ψ be finite. Then Assumption 2.1 holds for θ = |m(ψ 0 )| and θn = | mn (ψ 0 )| = 1 n n j=1 (X j -μ1,n (j)) e -iψ 0 j . Therefore Theorem 2.1 ensures that the subsampling based on Xn is consistent i.e., (i)-(iii) in Theorem 2.1 hold.

Testing if B(t, τ ) is periodic

In this section we discuss the problem of testing periodicity of the autocovariance function in the class of PC time series with known period T . Note that in such a class of time series we have µ 2 (t) ≡ 0. Moreover, μ1,n (t) is the estimator of the periodic function µ(t) and is given by formula (5) with Ψ = Ψ 1 = {2kπ/T : k = 0, 1, 2, . . . , T -1}. To simplify the notation we will write μn (t) instead of μ1,n (t). We take τ ∈ Z and we consider the following test

H 0 : B(t, τ ) is constant in t ∈ Z vs H 1 : B(t, τ
) is periodic in t with period T . Note that H 0 can be now equivalently rewritten as H 0 : |a(λ, τ )| = 0 f or all λ ∈ {1, 2, . . . , T -1}. Assume that τ ≥ 0. The estimator of a(λ, τ ) based on Xn has the form

ân (λ, τ ) = 1 n n-τ t=1 (X t -μn (t))(X t+τ -μn (t + τ ))e -iλt . (6)
Therefore we consider the following test statistics

U n = √ n λ∈A |â n (λ, τ )|,
where

A = {2kπ/T, k = 1, . . . , T -1}. Large values of U n suggest that H 1 is true.
To approximate quantiles of U n , we show subsampling consistency.

Theorem 3.2 Assume that conditions A2 and A3 are fulfilled. Then Assumption 2.1 holds for θ = λ∈A |a(λ, τ )| and θn = λ∈A |â n (λ, τ )|, where ân (λ, τ ) is given by (6). Therefore Theorem 2.1 ensures that the subsampling based on Xn is consistent i.e., (i)-(iii) in Theorem 2.1 hold.

A real data example

Below we present the real data application of the testing problem described in Section 3.2. For that purpose we use data published by Eurostat, concerning monthly industrial production (percentage change on previous period) for Mining and quarrying; manufacturing; electricity, gas, steam and air conditioning supply in three countries: France, Germany and Italy. In the first row of Figure 1 are visualized gross data i.e., data that were not modified in any way before publishing (neither seasonally nor calendar adjusted). The length of the analyzed sequence for each country is n = 288 (January 1992 -December 2015).

The considered data are known to contain two types of periodic structures. The first one represents the typical seasonal effect with period length equal to 12 months. The second one is a consequence of the so-called trading-day effect or the calendar effect, which is caused by a different number of working days in each month. Note that, each February in non-leap-years has the same number (equal to 20) of working days. For other months, this number varies from month to month. For the Georgian calendar this additional deterministic source of variability (well characterized for macroeconomic time series) is periodic with period equal to 400 years. But only a few frequencies corresponding to this periodic phenomenon are predominant (for more details see chapter 11.2 in [START_REF] Bell | Economic Time Series: Modeling and Seasonality[END_REF]). Thus, following the notation introduced in Section 3.2 we have that Ψ 1 ⊆ {2kπ/12, k = 0, . . . , 11}. Moreover, it is known that for monthly data the predominant trading-day frequency ψ 0 ∼ = 2.18733. We perform the test introduced in Section 3.2 to confirm that the frequency ψ 0 is significant i.e., that |m(ψ 0 )| = 0. Thus, we test the hypothesis:

ψ 0 ∈ Ψ 2 against ψ 0 ∈ Ψ 2 . Note that ψ 0 ∈ Ψ 2 is equivalent to |m(ψ 0 )| = 0.
In our approach we are using data obtained after removing the seasonal effect. One may note that in the case of the raw data, the ψ 0 is very close to the frequency ψ = 2π/3 ∈ Ψ 1 . To visualize this fact in the second and in the third row of Figure 1 we present the estimated values of T n = √ n| mn (ψ)|, ψ ∈ [0, π] obtained using the raw data X n and partially demeaned data Xn , respectively. The observed effect could have a serious impact on the results of our test if we used simply the raw data. To check that we decided to apply our methodology not only to the data obtained after removing the seasonal effect Xn , but also to the raw data X n (based on [Lenart, 2013]). Below we discuss results of both approaches. We set the significance level α = 5% and subsample length b ∈ {15, 20, . . . , 40}. The obtained values of the test statistics for the raw data are: 0.397 (Germany), 0.371 (France) and 0.497 (Italy). In the case of partially demeaned data the corresponding values are a bit lower: 0.382 (Germany), 0.303 (France) and 0.411 (Italy). For each considered country values of the subsampling quantiles are the highest for b = 15 and decrease when b increases. For example in the case of Germany taking b = 15 we got a subsampling quantile equal to 0.117 (raw data) and 0.06 (partially demeaned data). Corresponding values for b = 40 are: 0.094 and 0.042. Moreover, the raw data quantiles have much higher values than those for the partially demeaned data. As a result in the case of Italy values of the subsampling quantiles for the raw data are close to the values of the test statistics. Additionally, for Italy when the raw data and b = 15 are considered, the null hypothesis is not rejected i.e., at a significance level of 5%, frequency ψ 0 is not detected (p-value 0.0511). For b ∈ {20, 25, 30, 35, 40} the corresponding p-values are equal to 0.026, 0.0114,0,0,0 and hence H 0 is rejected. Simultaneously for the partially demeaned data, independently of the chosen country and the subsample length b, the null hypothesis is strongly rejected (p-values are equal to 0), which means that at the considered significance level we always detect the frequency ψ 0 . Thus, using our subsampling approach allows for better frequency detection.

Appendix

Proof of Theorem 2.1. The proof is similar to the one proposed by [START_REF] Politis | Subsampling[END_REF]] and hence we skip most of the details, indicating only some major differences. Let q = n -b + 1 and

U n (x) = 1/q q t=1 1{τ b ( θ n,b,t - θ) ≤ x}. Notice that Ln,b (x) = 1/q q t=1 1{τ b [ θ n,b,t -θ] + τ b [ θn,b,t -θ n,b,t ] - τ b [ θn -θ] ≤ x}.
We show that for any x ∈ R and > 0

U n (x -)1{E n } ≤ Ln,b (x)1{E n } ≤ U n (x + ), (7)
where

E n = max 1≤t≤q |S n,b,t | ≤ and S n,b,t = τ b [ θn,b,t -θ n,b,t ] -[ θn -θ] for 1 ≤ t ≤ q. If 1{E n } = 0, then (7) is obvious. If 1{E n } = 1 then max 1≤t≤q |S n,b,t | = τ b max 1≤t≤q [ θn,b,t -θ n,b,t ] -[ θn -θ] ≤ ,
which means that uniformly in 1 ≤ t ≤ q we have S n,b,t ∈ [-, ]. Thus, for any 1 ≤ t ≤ q we get that

1{τ b ( θ n,b,t -θ) ≤ x -} ≤ 1{τ b ( θ n,b,t -θ) ≤ x -S n,b,t } = 1{τ b [θ n,b,t -θ] + τ b [ θn,b,t -θ n,b,t ] -τ b [ θn -θ] ≤ x}, (8)
1{τ b ( θ n,b,t -θ) ≤ x + } ≥ 1{τ b ( θ n,b,t -θ) ≤ x -S n,b,t } = 1{τ b [θ n,b,t -θ] + τ b [ θn,b,t -θ n,b,t ] -τ b [ θn -θ] ≤ x}. (9)
Summing inequalities ( 8), ( 9) over t = 1, 2, . . . , q we get (7).

In the next step we show that

P (E n ) → 1. Note that max 1≤t≤q |S n,b,t | ≤ max 1≤t≤q τ b | θn,b,t -θ n,b,t | + τ b | θn -θ|.
Using assumptions (iii) and (i) we obtain that max Proof of Theorem 3.1. Notice that under A1 the condition (ii) of Assumption 2.1 holds by the same arguments as were used in the proof of Theorem 2.2 in [Lenart, 2013]. To show (i) note that

1 n n t=1 (X t -μ1,n (t))e -iψ 0 t -m(ψ 0 ) = 1 n n t=1 (X t -µ 1 (t))e -iψ 0 t -m(ψ 0 ) + 1 n n t=1 (µ 1 (t) -E[μ 1,n (t)])e -iψ 0 t + 1 n n t=1 (E[μ 1,n (t)] -μ1,n (t))e -iψ 0 t . ( 10
)
The second and the third summand on the right-hand side we denote by p n and q n , respectively. We have

√ np n = 1 √ n n t=1 ψ∈Ψ 1   m(ψ) -   1 n n j=1 ρ∈Ψ m(ρ)e iρj e -iψj     e iψt e -iψ 0 t = 1 √ n ψ∈Ψ 1   1 n ρ∈Ψ\{ψ} m(ρ) e i(ρ-ψ) (e in(ρ-ψ) -1) e i(ρ-ψ) -1   e i(ψ-ψ 0 ) (e in(ψ-ψ 0 ) -1) e i(ψ-ψ 0 ) -1 . Therefore | √ np n | = O(n -3/2 ). Moreover, note that | √ nq n | is equal to 1 √ n   ψ∈Ψ 1   1 n n j=1 µ(j)e -iψj   - ψ∈Ψ 1   1 n n j=1 X j e -iψj     e i(ψ-ψ 0 ) (e in(ψ-ψ 0 ) -1) e i(ψ-ψ 0 ) -1 ≤ 1 √ n ψ∈Ψ 1   1 n n j=1 |µ(j) -X j |   2 |e i(ψ-ψ 0 ) -1| .
Since E|µ(j) -X j | is bounded by some constant independent of n uniformly at j ∈ Z, we get that E| √ nq n | = O(n -1/2 ) and finally that √ nq n p -→ 0. To show (iii) we use the inequality ||x| -|y|| ≤ |x -y| (which holds for any x, y ∈ C). Then for any 1

≤ t ≤ n -b + 1 the term τ b | θn,b,t -θ n,b,t | can be bounded from above by 1 √ b t+b-1 j=t (X j -μ1,n (j))e -iψ 0 j - t+b-1 j=t (X j -µ 1 (j))e -iψ 0 j = 1 √ b t+b-1 j=t   ψ∈Ψ 1 m(ψ)e iψj - ψ∈Ψ 1 1 n n k=1 X k e -iψk e iψj   e -iψ 0 j ≤ 1 √ b ψ∈Ψ 1 m(ψ) - 1 n n k=1 X k e -iψk 2 |e i(ψ-ψ 0 ) -1| . Since |m(ψ) -1/n n k=1
X k e -iψk | does not depend on t and tends to 0 in probability as n -→ ∞ for each ψ ∈ Ψ 1 , we get that (iii) holds.

Proof of Theorem 3.2. Notice that under A2 and A3 the condition (ii) of Assumption 2.1 holds by the same arguments as were used in the proof of Theorems 1 and 2 in [START_REF] Lenart | Subsampling in testing autocovariance for periodically correlated time series[END_REF]. To show (i) one needs to decompose √ n( θn -θ) similarly to (10) and then to show the convergence to 0 in probability of the corresponding terms. For that purpose let us denote

p n = 1/ √ n n-τ t=1 (X t -µ(t))(µ(t + τ ) -μn (t + τ ))e -iλt , q n = 1/ √ n n-τ t=1 (µ(t)-μn (t))(X t+τ -µ(t+τ ))e -iλt and r n = 1/ √ n n-τ t=1 (µ(t)- μn (t))(µ(t + τ ) -μn (t + τ ))e -iλt .
Note that using the same arguments as in the proof of Theorem 2.1 in [Lenart, 2013] (under the same assumptions) and results of [Kim, 1994] we get that mn (ψ) -m(ψ) 2 = O(n -1/2 ). ( 11) Thus,

E|p n | ≤ ψ∈Ψ m(ψ) -mn (ψ) 2 1 √ n n-τ t=1 (X t -µ(t))e -iλt+iψ(t+τ ) 2 = o(1).
Similarly E|q n | -→ 0 as n -→ ∞ and hence p n , q n -→ 0 in probability as n -→ ∞. Finally, E|r n | is bounded by

ψ 1 ∈Ψ ψ 2 ∈Ψ m(ψ 1 ) -mn (ψ 1 ) 2 m(ψ 2 ) -mn (ψ 2 ) 2 1 √ n n-τ t=1 e -iλt+iψ 1 t+iψ 2 (t+τ ) = ψ 1 ∈Ψ ψ 2 ∈Ψ m(ψ 1 ) -mn (ψ 1 ) 2 m(ψ 2 ) -mn (ψ 2 ) 2 O( √ n) = o(1).
To get (iii) note that using similar arguments as in the proof of Theorem To finish the proof of the theorem it is enough to show that the each term of the sum on the right-hand side of ( 12) converges to 0 in probability as n -→ ∞. For that purpose let us take any ψ 1 , ψ 2 ∈ Ψ. If |m(ψ 1 )m(ψ 2 )| = 0 (without loss of generality we assume that |m(ψ 1 )| = 0), then O( b/n) √ n| mn (ψ 1 )|| mn (ψ 2 )| -→ 0 in probability as n -→ ∞. The last convergence follows from convergence of √ n mn (ψ 1 ) in distribution to Gaussian random variable (see Theorem 2.1 in [Lenart, 2013]), convergence O( b/n)| mn (ψ 2 )| p -→ 0 as n -→ ∞, the continuous mapping theorem and Slutsky's lemma. If |m(ψ 1 )m(ψ 2 )| = 0, then from Theorem 2.1 in [START_REF] Lenart | Seasonality Revisited -statistical Testing for Almost Periodically Correlated Stochastic Processes[END_REF], the multivariate delta method and Slutsky's lemma we have that O( b/n) √ n| mn (ψ 1 ) mn (ψ 2 ) -m(ψ 1 )m(ψ 2 )| p -→ 0 as n -→ ∞.

Figure 1 :

 1 Figure 1: First row: Monthly industrial production (percentage change on previous period, raw data -neither seasonally nor calendar adjusted, source: Eurostat) from Jan. 1992 to Dec. 2015. Second and third rows: values of test statistics T n obtained respectively for the raw data X n and for data got after subtracting estimates of seasonal means Xn .

  we get that P (E n ) → 1.

ψ 1

 1 3.1, expression τ b | θn,b,tθ n,b,t | can be bounded from above μ n (j + τ ) -µ(j + τ )]e -iλj + (j)μ n (j + τ ) -µ(j)µ(j + τ )]e -iλj , for any 1 ≤ t ≤ n -b + 1. Let us denote summands above by a n , b n and c n , respectively and q = n -b + 1. Note that max1≤t≤q τ b | θn,b,t b -θ n,b,t | ≤ max ψ∈Ψ | mn (ψ) -m(ψ)|M b, a.s.,where the last inequality follows from assumption A2. Hence by (11)max 1≤t≤q a n 2 ≤ ψ∈Ψ mn (ψ) -m(ψ) 2 O(b 1/2 ) = O( b/n) → 0.This means that max 1≤t≤q a n -→ 0 in probability as n -→ ∞. The corresponding convergence for b n can be obtained analogically. Finally, ∈Ψ ψ 2 ∈Ψ | mn (ψ 1 ) mn (ψ 2 ) -m(ψ 1 )m(ψ 2 )|. (12)

Acknowledgments

The authors would like to express their sincere gratitude to the anonymous referees for their valuable suggestions. Lukasz Lenart was supported by the Polish National Science Center based on decision number DEC-2013/09/B/HS4/01945. Anna Dudek has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie Grant Agreement No. 655394.

Subsampling for nonstationary time series with non-zero mean function