Zhen Qu

Dianyin Hu

Rongqiao Wang

Investigation on Dynamic Meshing Process and Factors Influencing Root Crack Propagation Trajectory of a Representative Aero-engine Gear Pair

Keywords:

The gear in an aero-engine experiences alternate loads therefore the stresses on gear root change significantly in meshing process. As a consequence, fatigue crack initiates on the gear root. This paper examined meshing process and crack propagation trajectory of the aero-engine gear based on explicit dynamics and linear elastic fracture mechanics. The crack on gear root affects meshing process and causes additional noise and vibration.The meshing impact and the peak of stress distribution of the gear pair is also magnified. Furthermore, the root crack propagation can cause the rim fracture and tooth fracture of gear system. The occurrence of rim fracture increases as the backup ratio (i.e., rim thickness divided by tooth height) decreases and also increases as the initial crack location is moved down the root of the tooth; the orientation of initial crack has little influence on crack propagation trajectory.

INTRODUCTION

Gear transmission, which is widely used to transfer the torque and force between any two or more axles, is the most important form in mechancial drive. Main accessory systems of the engine, such as fuel system, lbrication system and hydraulic system, are all powered by the engine rotor though the gear transmission system. The tooth of the gear endures the fatigue loads(meshing force), which causes the fatigue crack on the root of the tooth during the meshing process.

Numerous studies on dynamic meshing process and the crack propagation trajectory of the gear have been carried out.

Hu analyzed the dynamic characteristic (time-varying meshing stiffness, acceleration, error transfer, mesh force, friction) of gear system with root crack using the lumped parameter dynamic model. [START_REF] Hu | Fault Detection of Gears with Different Root Crack Size Using Wavelet[END_REF][START_REF] Hu | Influence of Different Gear Generation Methods on Contact Stress Analysis of Spur Bevel Gears in PSD[END_REF][START_REF] Wang | The Mathematical Model of Spiral Bevel Gears-A Review[J][END_REF] Wu S used the statistical methods to study the vibration response of the gear pair with simplified linear crack. [START_REF] Wu | Simulation of Spur Gear Dynamics and Estimation of Fault Growth[END_REF] Although the dynamic meshing process of the crack gear is analyzed, the crack propagation path is not further considered.

Lewicki from fracture mechanics group of Cornell studied the influencing factors (initial position of the crack, initial angle of the crack, geometic construction of the gear and the rotational speed of the gear) of the crack propagation trajectory. [START_REF] Lewicki | Gear Crack Propagation Path Studies[END_REF][START_REF] Lewicki | Effect of Speed (Centrifugal Load) on Gear Crack Propagation Direction[END_REF][START_REF] Lewicki | Rim Thickness Effects on Gear Crack Propagation Life [J][END_REF][START_REF] Lewicki | Crack propagation studies to determine benign or catastrophic failure modes for aerospace thin-rim gears[END_REF] Liu simulated the trend of crack propagation of the maximum stress area of the gear tooth using ABAQUS. [START_REF] Shuang | Simulation on propagation of Involute Gear Tooth Based on ABAQUS[END_REF] Curà F calculate the impact of the initial position and angle of crack and the ratio of hight and the thickness of tooth by the method of extend finite element method. [START_REF] Curà | Investigation about crack propagation paths in thin rim gears[END_REF]. Spievak L E investigated the crack trajectory prediction and fatigue life for a spiral bevel pinion using Finite Element Method. [START_REF] Spievak | Simulating fatigue crack growth in spiral bevel gears[END_REF] Based on the aforementioned analysis, this paper focuses on a respective aero-engine gear pair in the accessory gearbox and makes the comprehensive consideration of dynamic meshing process and crack propagation trajectory.

METHODS

In the study of the dynamic meshing process of the gear pair with root crack, based on the method of explicit dynamics analysis, the time-varying parameters(the rotational speed, engaging force and Von Mises stress) of the cracked gear pair are calculated and campared with the healthy one. And the analyzing and processing these time-varying signals provides the basis to the diagnosis of the gear pair. If there appears initial crack on the root of gear tooth, it will grow during meshing process. The crack propagation trajectory differs on different initial condition and different gear design. So, using the method of finite element and the code of franc2d [START_REF] Wawrzynek | FRANC2D: A two-dimensional crack propagation simulator[END_REF], the influencing factors (initial position, angle of the crack and backup ratio) of crack propagation trajectory and corresponding different fracture forms(tooth fracture and rim fracture) are carried out. The entire process of crack propagation is divided into several steps, and the angle (θ ) of each crack propagation can be calculated by the following procedure. Firstly, with the method of J-integral [START_REF] Cherepanov | The propagation of cracks in a continuous medium[END_REF] [START_REF] Rice | A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks[END_REF] in the finite element analysis, the stress intensity factor and the stress distribution can be calculated using the code of franc2d.

[ dy T ds]

u JW x Γ ∂ = -∂ ∫  (1) 0 11 22 ij ij ij ij ij Wd ε σ ε σε = = ∫ (2)
Equation( 1) is the J-integral around the crack tip and W in equation( 2) is energy density which can be got by the product of stress and strain.T is the force acting on the integral boundary.

K JG E = = (3)
Because the result of J-integral is equal to the enegry release rate G [START_REF] Rice | Fracture Mechanics[END_REF] in linear facture mechanics, stress intensity factor K can be got according equation [START_REF] Wang | The Mathematical Model of Spiral Bevel Gears-A Review[J][END_REF].E is elasticity modulus of material.

Figure1. The crack propagation trajectory of polylines

Then, the vertical direction i θ of the maximum circumferential stress around crack tip is regared as the direction of crack propagation for the step i (i=1,2,…).(Figure1.) Thirdly, the crack extends a tiny distence along i θ and the first polyline can be achieved. The codes continue to calculate the direction( 1 i θ + ) for the next step and the second polyline. By cycle computing, finally, the crack propagation trajectory can be obtined which is composed of these tiny polylines.

The following study of crack propagation trajectory prediction of the representative gear provides recommended design principles in order to prevent the occurrence of fracture failure of gears. In short, they can serve as the guidelines when conducting the failure prediction and fault diagnosis of gear system.

The geometry model and mesh generation 3.

The geometric parameters of gear pair are listed as follows.(Chart1 20°T he gear pair studied in this paper is rotated around z axis in the x-y plane. In order to save calculation time and improve the computational efficiency, the 2d gear model which is the projection on the x-y plane is used in the calculation. In analysis of dynamic meshing process and crack propagation trajectory, the 2d gear pair model(Figure2) and the 2d partial gear model(Figure3) with three teeth is used separately. In the FEM analysis, the elasticity modulus is 194GPa, the poission ratio is 0.297 and the density is 7.86e-9t/mm 3 .Besides, the element type is plane stress element with midside node. (driving gear) In the process of generation of elements of the gear pair, each gear is divided into two parts. So the whole finite element model is divided into four parts shown in Figure2, driving gear part (red element part), driving gear rigid part (green element part), driven gear part (blue element part) and driven gear rigid part (yellow element part). The rigid element parts of gear pair which can not be denformed are only used to transfer velocity to the driving part and resistance torque to driven part. The crack on the tooth root (Figure4.) is generated by dividing one node on the root of gear tooth into two nodes and its length depends on the side length of the element on the driving gear.

Analysis of dynamic meshing process 3.1 Contact dynamics method

During the process of gear meshing which involves impact and contact between gear teeth, the magnitude and direction of loads(engaging force) acting on the tooth change as time. When solving the contact problems, the penalty method is used for its simple principle, easy programing, low possibility of causing hourglass effect, less numerical noise and the high precision in the method of momentum conservation. [START_REF] Bai | LS DYNA3D theory foundation and example analysis[END_REF] fk δ = (4)

2 s f KA k V = ( 5
)
The penalty method assumes that there is linear spring between the contact points if they have penetrated each other and the spring force is expressed as Equation( 4), where f is the spring force, k is the contact stiffness. The contact stiffness as Equation [START_REF] Lewicki | Gear Crack Propagation Path Studies[END_REF], where s f is penalty factor(default value is 0.1), K is the bulk modulus of contact elements, A is elements contact dimension and V is elements contact volume.

The aim of dynamic analysis of meshing process is to compute the mechanical and time parameters of each time-step The dynamic parameter of next time-step is related to the meshing state of pervious time-step. So on the basis of the previous meshing state, an incremental solution procedure can be used to solve the problem. First, we should determine the time step t ∆ in accordance with the reference [START_REF] Wu | Contact Finite Element Analysis and Application Study on Dynamic Meshing Characteristics of Gear Transmission System[END_REF] and the whole analysis of meshing process can be dispersed into several time series(0, ,

+ ++ +             = --                      ; (7
+ ++ +             = --                      ; ( 8
)
tt t tt t g g gp p p u u uu u u ++ = -= -  ; ( 9
)
Subscripts g and p are for driving gear and driven gear. and t u (t=0) can be determined by the initial condition of the gear engagement. According the penalty method, the nodal contact force vector can be written as ( )

 = -- - +   + =++      ( 10
)
where C N is interpolation function, µ is coefficient of sliding friction, T u is tangential displacement incerment,
( )

1 2 tt t t tt u u u ut +∆ + = ++ ∆      (12) 1 2 1 2 tt t tt tt tt LC tt F F F u ut u tC +∆ +∆ +∆ +∆ + - -+ ∆ = +∆ [] [] []( ) [ C ] [M] [ ]    (13) [] [] = tt tt F F Ku + +  [] (14)
Through equation ( 11) and ( 13), tt u +∆ and tt u +∆  can be figured up and then tt u +∆  can be calculated subsenquently by Equation [START_REF] Wawrzynek | FRANC2D: A two-dimensional crack propagation simulator[END_REF]. Then the internal force, contact force and nodal foce of next increment time step can be calculated.

The analysis of the meshing calculation result

The loading curves are shown in Figure5(a) and (b).The rotating speed(1608.6rad/s) is applied on the driving gear rigid part and the burden torque (141900N•mm) is applied on the driven gear rigid part. The value of rotating speed and torque match actual parameters of a genuine engine. Automatic surface-to-surface contact is selected as the contact principle.

Analysis of the rotating speed

The rotating speed applied on driving gear rigid part starts at The trends of the output rotating speed of both the healthy gear and the cracked one coincide with the input rotating speed of the driving gear(Figure5). After 0.01 second the rotating speed of the driving gear reaches the setting speed linearly and the rotating speed of the driven gear fluctuates around the theoretical rotating speed. By comparing the curve in Figure6 and Figure7, we can see that the fluctuation amplitude of the rotating speed of the healthy driven gear is almost stable after 0.04 second and the cracked gear is 0.065 second. The maximum rotation speed fluctuation amplitude of the cracked gear is about 750rad/s (point A in Figure7)which is almost ten times larger than the healthy one (75 rad/s) during meshing process. So, it can be seen that the vibration shock of the cracked gear pair is significantly enhanced in the process of engaging. Affected by the increase of vibration amplitude, the gear pair inevitably produces the abnormal noise during the meshing process, which has a seriously adverse impact on the normal stable transmission of gear pair.

The time domain data of 0.08 second to 0.1 second in Figure6 and Figure7 which have been stable are subjected to Fourier transform to obtain the curve in frequency domain shown in Figure8(a) and (b) which is calculated by LS-DYNA. The first peak on the left of the curve in Figure8(a) and (b) corresponds to the rotational frequency which is 298Hz(1608.8/2π) theoretically. The peaks of the middle part of the curve corresponds to the meshing frequency which is 7936.67Hz in theory and the possible natural frequency of the gear. The rightmost peak corresponds to the second frequency of the meshing frequency. Because of the crack on the root of the gear the proportion of the meshing frequency and the vibration amplitude is significantly changed. The crack on the root of the gear pair that causes damage to the structure of the gear pair intensifies the impact in the meshing process, which causes abnormal noise and vibration and greatly shorten the service life of the gear pair

Analysis of the engaging force and the element contact stress

The output of contact force and the node force come from the result of the LS-DYNA calculation, the engaging force is extracted shown in Figure9.

In the initial stage (0.01s-0.02s) of high-speed meshing, the curve of the meshing force fluctuates greatly with time, and it is obvious that there is an impact. And with the time goes by, the fluctuation amplitude of the engaging force decreases gradually. After about 0.04 seconds in Figure9(a) and 0.065 seconds Figure9(b), the mesh enters the stable transmission phase, which coincides with the result of the rotating speed in Figure8. The engaging force of the healthy driving gear is constant positive, which means that there is only unilateral meshing which means that only one side of the tooth is in the contact situation during the process but the engaging force of the cracked gear has both positive and zero value, based on which the engaging force of the driving gear should be decomposed into the xdirection(parallel to the direction of line of centers) and ydirection (perpendicular to the direction of line of centers) shown in Figure10.

(a) Healthy Driving Gear (b) Cracked Driving Gear Figure9. Engaging Force Compared with the meshing process of the healthy gear pair, the moment when the cracked gear pair enters the stable speed fluctuation which means that the vibration amplitude of speed is almost constant is obviously pushed backward. Moreover, the maximum peak-to-peak value of the fluctuation of engaging force of the cracked gear is 25KN, which is much larger than that of 5KN of healthy gear pair. The x-direction of engaging force shown in Figure10(a) exists both positive value and zero value indicates that there is a phenomenon of tooth off which means that there is a moment when the tooth of the gear is not in the contact condition in meshing process. The x-direction of the engaging force shown in Figure10(b) exists positive, zero and negative value indicates that there are not only unilateral and 'tooth off' meshing but also 'bilateral impact' meshing which means that both the two sides of one tooth have an impact on the adjacent sides of other tooth. The cause of the occurrence of tooth off and bilateral impact is that the ability of gear teeth to resist bending deformation is reduced and the bending stiffness of the gears is also reduced after the initiation of crack during the meshing process. The equivalent stress(Von Mises stress) of four chosen elements on the contact surface of gear tooth changing over time is shown in Figure11 during the meshing process of one tooth. By comparing (a) and (b) in Figure11, obviously, there are two peaks on the curve of the cracked gear while only one on the curve of healthy gear, which indicates that ,because of the crack on the root of the gear intensifying the vibration shock, bilateral impact occurs on the cracked gear pair while there is only unilateral meshing on the healthy gear pair. The maximum equivalent contact stress (Von Mises stress) of cracked gear is 660MPa, which is obviously larger than that of the healthy gear and this is also caused by the increased impact of cracked gear.

Analysis of the crack propagation trajectory

The initial crack on the root of gear will propagate during the meshing process. At the initial propagation stage, the crack on the gear is so tiny that could not cause the catastrophic failure of gear system. Acceleration sensor can be used to monitor data of the gear system and determine whether there is a failure in the system according the analysis of previous chapter. However, if the initial crack is not monitored at the initial propagation stage and the gear system continues to work, the initial crack will continue to propagate in different directions and cause the rim fracture or tooth fracture. The rim fracture leads to more serious failure in the gear system than tooth fracture, so it is necessary to analyze the possible crack propagation trajectory.

4.1Determination of load and boundary conditions

The literature [START_REF] Lewicki | Consideration of Moving Tooth Load in Gear Crack Propagation Predictions[END_REF] combined the numerical simulation and experimental verification to draw the conclusion that in the simulation of the tooth root crack propagation trajectory, the crack propagation path can be obtained sufficiently accurately with a static load applied to the gear at the position of the highest meshing point instead of the moving dynamic load. The value of the equivalent static load is determined according to the actual conditions of the gear pair. of the gear tooth is obtained by equation ( 7), ( 8) and ( 9).

* (z h ) m

a d = +  ( 16
)
2000 input t T F d = (17) 20 tan rt FF =   ( 18
)
The force is applied on the point A and the fix constraint is applied to the lower, left and right sides of the rim of the partial gear model.(Figure12)Experiments on the finite element model of the whole gear and partial gear pair with three tooth have been done ,which justifies that the constraints used on the partial gear pair has the same result with whole gear.

Figure12. Load and Boundary Condition

Influencing factor of crack propagation trajectory

The angle of crack propagation is obtained for each sub-step by calculating the stress intensity factor of crack tip.The local meshes of the crack region can be redivided by the 'element addition and deletion' in the FRANC program and the program doesn't sotp until the crack propagation trajectory reaches the boundary of the finite element model. Finally, the crack propagation trajectory can be obtained shown in Figure13. The simulation results of the crack propagation trajectory coincide with the simulated results in the literature [START_REF] Hu | Influence of Different Gear Generation Methods on Contact Stress Analysis of Spur Bevel Gears in PSD[END_REF] and the experimental results in literature [START_REF] Lewicki | Consideration of Moving Tooth Load in Gear Crack Propagation Predictions[END_REF] show the accuracy of the numerical simulation of the FRANC program in literature [START_REF] Hu | Influence of Different Gear Generation Methods on Contact Stress Analysis of Spur Bevel Gears in PSD[END_REF].

Figure13. Crack Propagation Trajectory

Influence of initial crack location on crack propagation trajectory

The initial crack is inserted at the root of the gear model from point 1 to point 8 shown in Figure14(a). The initial crack length is set to 0.3mm along the direction of the element. The crack is formed using the same method with finite element model in the dynamic calculation, but the length is shorter because the element size of the partial model with three teeth is smaller in order to obtain preciser crack propagation trajectory. The simulation of 8 crack propagation trajectories which is integrated into one figure is completed using the FRANC program 8 times shown in Figure14(b). The length of the crack trajectory in the rim region is gradually increased as the initial crack position moves down, and all the fracture form is the tooth fracture.

Influence of initial crack orientation on crack propagation trajectory

In the process of gear processing, there are subtle scratches and the concentration of notch stress on the surface of the material and the corrosion of the working environment. When the local stress concentration exceeds the breaking strength, the atomic bond breaks and causes the initial fatigue crack in different directions on the surface of the material. The difference of the initial direction of the crack has a certain impact on the crack propagation trajectory.

The gears are subjected to alternating loads during the working process and the root is the position of stress concentration where there are the dangerous sections. The position of the dangerous section is determined by many kinds of methods.In this paper, the Hofer 30 ° Hazard section method [START_REF] Xue | The calculation of the dangerous cross section of gear teeth based on conjugate surface principle[END_REF] is used. According to the results of a large number of gears opto-elastic experiments, the tangent points (point A and point B in Figure 16) are obtained by drawing the line tangent to the tooth profile at a 30 ° angle with the centerline. The position of connection line between point A and point B is regarded as the position of the dangous section of the gear. is changed(Figure18), the crack propagation trajectories is basically the same as those in Figure17. Crack propagation causes the tooth fracture, so the influence of the initial crack oriention on the crack propagation path can be used as the secondary design reference factors.

Influence of backup ratio on crack propagation trajectory

From the results of previous chapter, it can be deduced that when the rim thickness is thinned, the crack propagation may lead to rim fracture. When the backup ratio ( b m ) is 0.5 and 0.2, the crack propagation trajectory can be calculated with the same method and the result is shown in Figure19. With the thinning of the rim thickness, the form of tooth fracture caused by crack propagation is more likely to occur at rim, which causes a more serious catastrophic acciedent in the entire gear system compared with tooth fracture. For gears in the gearbox of aeroengine accessory gears, gear failure can lead to failure of the hydraulic system, the oil system, and the electronic control unit, which may cause the consequences of the crash. Therefore, in the process of design and strucutral optimization of the gears, the backup ratio should first be considered to ensure the safe and reliable operation of the gear system.

5.Conclusion

Based on the finite element program, this paper completes the dynamics simulation of the gear meshing and the prediction of crack propagation trajectory in the actual working condition, which provides the basis for the fault prediction and structural design of the gears.

The following conclusions are obtained from the study of this paper.

(1)The speed of the cracked gear fluctuates greater than that of the healthy one and is not easy to enter stable meshing stage.The crack changes the frequency domain characteristic of the rotating speed, which produces additional abnormal vibration and noise in meshing process.

(2)The crack on the root of the gear affect the stiffness of the gear, which can cause additional bilateral impact and disengagement. This leads to an increase in the instantaneous fluctuation of the engaging force, which can cause additional damage to the gear system.

(3) The position of the initial crack has a great influence on the crack propagation trajectory of the gear. The possibility of occurrence of rim fracture increases as the initial crack position moves down along the tooth root.

(4) The backup ratio of the gear structure also has a large effect on the crack propagation trajectory, and the possibility of occurrence of rim fracture is significantly increased with the decrease of backup ratio.

(5) The direction of the initial crack has little effect on the crack propagation trajectory, which can be used as a secondary reference factor in the design and optimization of the gear to ensure the safe and reliable operation of the gear system.

Figure2. 2d

 2d Figure2. 2d Model of Gear Pair

Figure3.

  Figure3. 2d Partial Model(driving gear)

  Matrix [M] , [] K and [C] are for mass matrix, stiffness matrix and dumping matrix. u  , u  and u are for the acceleration, velocity and displacement of the element node. The superscripts L and C of [] F are for the external nodal load and nodal contact force. t g K []and t p K []are respectively stiffenss matrix of the driving gear and driven gear at the reference mement which is determined by the position and configruation of gear pair. t p F and t g F are the internal force vector of driving gear and driven gear at t . p u and g u are incermental nodal displacement of drving and driven gear. t u 

  Loading CurveThe result of rotating speed of the driven gear, the engaging force and the element contact stress changing over time is analyzed in the following characters.

  Figure7. Rotating Speed of cracked gear pair

  GearFigure8. Amplitude-versus-frequency Curve

Figure10.

  Figure10. Decomposition of the Engaging Force of the Cracked Gear

Figure11.

  Figure11. Equivalent Stress of Elemen

F

  The input torque (T) of the driving gear is 184.9N•m by equation[START_REF] Lewicki | Effect of Speed (Centrifugal Load) on Gear Crack Propagation Direction[END_REF] , where the rotating speed (N) is 15361.3 r/min and input power (P) of the driving gear is 297.5kw. The transmission ratio (i) is 31/40 and the transmission efficiency is 0.99. The input torque input T can be calculated. ) and nominal radial force ( r F )

  (a)Initial Crack Location (b)Crack Propagation Trajectory Figure14. Influence of the Initial Crack Loaction on Crack Propagation Trajectory

Figure16.Figure18.

  Figure16. Position of the Dangerous Section

Figure19.

  Figure19. Influence of Backup Ratio on Crack PropagationTrajectory The crack propagation results of point 8 in Figure19(a) and point 5,6,7,8 in Figure19(b) are rim fracture, which verify the above inference. With the thinning of the rim thickness, the form of tooth fracture caused by crack propagation is more likely to occur at rim, which causes a more serious catastrophic acciedent in the entire gear system compared with tooth fracture. For gears in the gearbox of aeroengine accessory gears, gear failure can lead to failure of the hydraulic system, the oil system, and the electronic control unit, which may cause the consequences of the crash. Therefore, in the process of design and strucutral optimization of the gears, the backup ratio should first be considered to ensure the safe and reliable operation of the gear system.

  ⋅⋅⋅ + ∆ ⋅⋅⋅ ).Second, time incremental algorithm is induced in the motion equation set on the concept of virtual displacement principle.(Equation6-9)

	, tT ∆ ∆ tt t t t gg tt L tt g tt pp tt L tt p M u F M u F +∆ +∆ +∆ +∆ +∆ +∆  +  [ ] [] gg tt C t gg tt pp C t pp C u FF C u FF +∆ +∆ += t g g t p p Ku Ku  []   +- [] [] []  + +=  +-  [ ] [] []    [ ] [] [] (6)
			tt	00 11 11 t xx xx
	tt	u	00 11 11 = t t gg tt yy u yy
			tt	00 11 11 t zz zz

  )

			tt	00 22 22 t xx xx
	tt	u	00 22 22 = tt pp tt yy u yy
			tt	00 22 22 t zz zz

ACKNOWLEDGMENTS

The present work is the result of all joint efforts. I would like to express my deep gratitude to Professor Dianyin Hu and Professor Jianjun Wang.who helped me modify the paper. Also, I would like to experss my sincere gratitude to all the classmates who have discuessed with me. Their instructions and suggestions have helped broden my horizon.

Last but not least, I would like to experss my speial thanks to my parents, whose care and support motivate me to move on.


F Meng
email: meng@ph1.uni-koeln.de

Á Sánchez-Monge

P Schilke

M Padovani

A Marcowith

A Ginsburg

A Schmiedeke

A Schwörer

C Depree

V S Veena

Th Möller

The physical and chemical structure of Sagittarius B2 V. Non-thermal emission in the envelope of Sgr B2

Keywords: Stars, formation -Stars, massive -Radio continuum, ISM -Radio lines, ISM -ISM, clouds -ISM, individual objects, Sgr B2 Previous studies

Context. The giant molecular cloud Sagittarius B2 (hereafter Sgr B2) is the most massive region with ongoing high-mass star formation in the Galaxy. In the southern region of the 40-pc large envelope of Sgr B2, we encounter the Sgr B2(DS) region, which hosts more than 60 high-mass protostellar cores distributed in an arc shape around an extended Hii region. Hints of non-thermal emission have been found in the Hii region associated with Sgr B2(DS). Aims. We seek to characterize the spatial structure and the spectral energy distribution of the radio continuum emission in Sgr B2(DS). We aim to disentangle the contribution from the thermal and non-thermal radiation, as well as to study the origin of the non-thermal radiation. Methods. We used the Very Large Array in its CnB and D configurations, and in the frequency bands C (4-8 GHz) and X (8-12 GHz) to observe the whole Sgr B2 complex. Continuum and radio recombination line maps are obtained. Results. We detect radio continuum emission in Sgr B2(DS) in a bubble-shaped structure. From 4 to 12 GHz, we derive a spectral index between -1.2 and -0.4, indicating the presence of non-thermal emission. We decomposed the contribution from thermal and non-thermal emission, and find that the thermal component is clumpy and more concentrated, while the non-thermal component is more extended and diffuse. The radio recombination lines in the region are found to be not in local thermodynamic equilibrium but stimulated by the non-thermal emission. Conclusions. Sgr B2(DS) shows a mixture of thermal and non-thermal emission at radio wavelengths. The thermal free-free emission is likely tracing an Hii region ionized by an O7 star, while the non-thermal emission can be generated by relativistic electrons created through first-order Fermi acceleration. We have developed a simple model of the Sgr B2(DS) region and found that first-order Fermi acceleration can reproduce the observed flux density and spectral index.

Introduction

The giant molecular cloud Sagittarius B2 (hereafter Sgr B2) is the most massive (∼ 10 7 M , see e.g., Goldsmith et al. 1990) region with ongoing high-mass star formation in the Galaxy, and has a higher density (> 10 5 cm -3 ) and dust temperature (∼50-70 K) compared to other star forming regions in the Galactic plane (see e.g., Ginsburg et al. 2016;Schmiedeke et al. 2016;Sánchez-Monge et al. 2017). Sgr B2 is located at a distance of 8.34 ± 0.16 kpc, at only 100 pc in projection from the Galactic center (Reid et al. 2014) 1 . In the central ∼ 2 pc, there are the two well-known hot cores Sgr B2(N) and Sgr B2(M) (see e.g., Schmiedeke et al. 2016;Sánchez-Monge et al. 2017), which con-tain at least 70 high-mass stars (see e.g., Gaume et al. 1995;De Pree et al. 1998, 2014). A larger envelope surrounding the two hot cores has a radius of 20 pc and contains more than 99% of the total mass of Sgr B2 (Schmiedeke et al. 2016, see also their Fig. 1 for a sketch of this region). The relatively high temperatures, densities, and pressure, as well as its proximity to the Galactic center make Sgr B2 a good target for studying highmass star formation in extreme environments. With this in mind, we have started a project to characterize the physical and chemical properties via observations and modeling (Schmiedeke et al. 2016;Sánchez-Monge et al. 2017;Pols et al. 2018;Schwörer et al. 2019). The current work is the fifth paper in the series. While the first observational studies focused on the two central hot cores Sgr B2(N) and Sgr B2(M), in this work we investigate the physical properties of the envelope they are embedded in.

regions Sgr B2(N) and (M). Martín-Pintado et al. (1999) report the presence of filament-, arc-and shell-shaped dense gas likely produced by stellar feedback. Recent observations with ALMA at 3 mm revealed more than 200 high-mass protostellar cores distributed throughout the envelope (Ginsburg et al. 2018). A remarkable feature is the presence of about 60 dense cores in the southern region, following an arc-shape distribution. This group of cores is spatially coincident with the arc-shaped feature detected in NH 3 by Martín-Pintado et al. (1999), which is referred to as Deep South (hereafter DS) by Ginsburg et al. (2018). At the position of Sgr B2(DS), Mehringer et al. (1993) report the detection of radio continuum and radio recombination line (RRL) emission (see also LaRosa et al. 2000;Law et al. 2008). Such continuum emission at radio wavelengths along with RRLs is usually related to thermal ionized gas from Hii regions (e.g., [START_REF] Kurtz | Star Workshop III: The Earliest Phases of Massive Star Birth[END_REF][START_REF] Kurtz | Massive Star Birth: A Crossroads of Astrophysics[END_REF]. However, along with the thermal emission detected in previous studies, non-thermal emission, has also been found in Sgr B2 due to relativistic electrons [START_REF] Larosa | [END_REF]Hollis et al. 2007;Protheroe et al. 2008;Jones et al. 2011). Yusef-Zadeh et al. (2007a, 2013, 2016) studied the presence of non-thermal emission in the Galactic Center region including Sgr B2. Particularly, observations at 255 MHz, 327 MHz, and 1.4 GHz (Yusef-Zadeh et al. 2007b) revealed non-thermal emission in Sgr B2(M). The thermal and non-thermal contribution to the radio continuum emission can be distinguished by their spectral energy distribution (hereafter SED). The SEDs of both thermal and non-thermal emission can be described by power-laws S ν ∝ ν α , where α is the so-called spectral index, which varies from -0.1 to +2 for thermal emission (see e.g., Sánchez-Monge et al. 2013) and becomes significantly negative for non-thermal emission (e.g., α = -0.8, see Platania et al. 1998).

In this paper, we present Very Large Array (VLA) observations of the whole Sgr B2 region in the frequency regime 4-12 GHz, with configurations BnC and D. Thus, this study focuses on spatial scales from 0.2-10 pc. The higher sensitivity and resolution compared to previous studies allow us to study in detail the presence and properties of radio-continuum sources, both thermal and non-thermal, throughout Sgr B2. Although the observations cover the whole envelope of Sgr B2, we pay special attention to Sgr B2(DS). In Sect. 2 we describe the observations as well as the data reduction process. Section 3 shows the results, and the study of the spectral index. In Sect. 4 we aim at decomposing the thermal and non-thermal components dominating the emission at radio wavelengths, while in Sect. 5 we discuss on the origin of the non-thermal emission. Finally, we summarize our study in Sect. 6.

Observations and data reduction

We used the VLA in its CnB and D configurations to observe the entire Sgr B2 complex in the frequency bands C (4-8 GHz) and X (8-12 GHz). The observations with the CnB configuration were conducted on May 3 and 5, 2016 (project 16A-195). On February 22 and 23, 2017, the D configuration was employed for the observations (project 17B-063). The continuum emission was observed by combining a total of 64 spectral windows with a bandwidth of 128 MHz each. Alongside the continuum observations, high-resolution spectral windows were used for line observations. Eighteen RRLs in the frequency range from 4 to 12 GHz were observed with a spectral resolution between 31.25 and 62.5 kHz (1-2 km s -1 ). In Table 1, we list the rest frequencies of the 18 RRLs. We used the mosaic mode to cover the whole extent of Sgr B2 (∼ 20 × 20 ). A total of 10 and 18 pointings, with primary beam sizes of 7.5 and 4.5 for the C and X Notes. (a) The frequencies of the stacked RRLs are labels corresponding to the average frequency of the stacked lines. These values do not correspond to actual transition frequencies, and are only used in the excitation analysis of Sect. 4.2.

bands, respectively, were used. The quasar 3C286 was used as a flux and bandpass calibrator. The SED of 3C286 from 0.5 to 50 GHz was measured by Perley & Butler (2013), with a flux of 5.059 ± 0.021 Jy at 8.435 GHz and a spectral index of -0.46.

The quasar J1820-2528, the flux of which is 1.3 Jy in the C and X bands, was used as phase calibrator. The calibration was done using the standard VLA pipelines provided by the NRAO2 . Calibration and imaging were done in Common Astronomy Software Applications (CASA) 4.7.2 [START_REF] Mcmullin | Astronomical Data Analysis Software and Systems XVI[END_REF]). The calibrated measurement sets of the CnB and D array data were concatenated to improve the uv sampling. Self-calibration was performed to reduce imaging artifacts and improve the final sensitivity of the image. Three loops of phase-only selfcalibration were conducted, with solution time intervals of two seconds, equal to the integration time of the observations. The self-calibrated measurement sets were used for imaging. All the pointings of the mosaic in each band were primary beam corrected and the mosaic was imaged using the CASA task tclean. With a robust factor of 0.5, the images of the C and X bands have synthesized beams of 2.7 ×2.5 , with a position angle (PA) of -82 • and 1.8 ×1.5 (PA = 76 • ), respectively. The PA is defined positive north to east. Under such resolutions, the root mean square (RMS) noise of the C and X band images are 0.4 mJy/beam and 0.2 mJy/beam respectively. This RMS noise is limited by dynamic range effects due to the bright emission in the region. We achieve a dynamic range of about 4000 and 2000 for the X-and C-band images, respectively. To investigate the SED over the whole 4-12 GHz range, the measurement sets of the C and X bands were divided into 12 frequency ranges. Thus, 12 images from 4 to 12 GHz were obtained. The uv coverage was restricted to 0.6-50 kλ for each frequency range to ensure that F. Meng et al.: Non-thermal emission in the envelope of Sgr B2 every image is sensitive to same spatial scales. The 12 images are convolved to a final circular Gaussian beam of 4 . In addition to the continuum emission, we also imaged the 18 RRLs. In order to increase the signal-to-noise ratio, we stacked the neighboring RRLs to produce a final set of four stacked RRLs (see Table 1 for details). All the RRL images were resampled to a common spectral resolution of 2.5 km s -1 . Finally, and with the aim of compensating for the low sensitivity (3 mJy/beam per 2.5-km/s channel), the images of the four stacked RRLs are smoothed to a resolution of 8 , resulting in a final sensitivity of 1 mJy/beam per 2.5-km/s channel.

Results

In this section we discuss the distribution of the radio continuum emission in the Sgr B2 region, characterize the spectral index, and study the RRL emission. We pay special attention to the Sgr B2(DS) region.

Ionized gas in Sgr B2(DS)

The images of the Sgr B2 region in the C and X bands are shown in Fig. 1. Some relevant objects such as the regions Sgr B2(N), (M) and (S) are indicated, as well as source V and the Hii region AA (see Mehringer et al. 1993). The DS region is the main focus of this article. It appears as a bubble-like structure with an outerdiameter of ∼1.5 pc. The thickness of the bubble edge is ∼0.3 pc, with the eastern part being stronger than the western region of the bubble, in both X and C-band images. Filaments and arcs are found along the edge of the bubble, while the emission towards the center is just at the RMS noise level.

The radio continuum emission in DS is spatially connected to the high-mass dust cores revealed by previous ALMA observations at millimeter wavelengths (Ginsburg et al. 2018). As shown in Fig. 2, more than 60 dust cores are located at the outskirts of the bubble-like structure, with a larger population of dense cores to the east.

Spectral index analysis

The C and X band data, from 4 to 12 GHz, are divided into 12 tomographic maps and convolved to a common angular resolution of 4 . To ensure that each of the 12 channel maps is sensitive to similar spatial scales, we have set the uv limit to the common range 0.6-50 kλ. Moreover, we have checked that large-scale background emission does not affect the spectral analysis. For this, we have inspected intensity profiles of the 12 channel maps, and found the background emission to be zero (see Appendix A and Fig. A.1). In Fig. 3, we show the portion of the map corresponding to Sgr B2(DS). The intensity decreases significantly from 4 GHz to 12 GHz. We analyzed the SED of the radio emission in DS by fitting a power law to the 12 tomographic maps. We define the variation of flux with frequency as S ν ∝ ν α , where S ν denotes the flux density, ν stands for frequency, and α is the spectral index. The power-law fitting is conducted pixel by pixel for the whole map. For each pixel, the spectral index is fit only if the emission in all the 12 tomographic maps is above 3σ, otherwise the pixel is masked. To avoid effects of possible artifacts, the fitting minimizes a loss function r(z) = √ 1 + z/0.01 -1, where z is the square of the difference between the fit power-law function and the observed intensity in the 12 frequency ranges. The uncertainty of α is obtained from the corresponding diagonal element of the covariance matrix of the fitting.

The fit spectral index map of the whole Sgr B2 region, as well as the map of uncertainty of α are shown in Fig. 4. The hot cores N, M, as well as other sources in the whole Sgr B2 region, except DS, show spectral indices between -0.01 and 2, which are consistent with thermal emission from Hii regions. The central parts in N and M have spectral index values greater than 1, which indicates that the emission is optically thick. At the edges of these hot cores we determine flatter SEDs (e.g., α = -0.1), suggesting that the centimeter emission at the edges of these Hii regions becomes optically thin.

In contrast to most regions in Sgr B2, the major part of DS shows negative spectral index values (see also Fig. 5). We have highlighted the SED of six selected positions in the DS region and shown them in Fig. 5. Source AA is included as reference for the spectral index, since it is a known Hii region with optically thin free-free emission (Mehringer et al. 1993). As displayed in Fig. 5, source AA has α = 0.01 ± 0.06, consistent with optically thin free-free emission. The other five spots have negative spectral indices, ranging from -0.4 to -1.2. Such values indicate the existence of non-thermal emission in DS. As shown in the map, the distribution of the non-thermal emission extends to more than 1 pc and forms arc structures.

With the goal of confirming the spectral index values derived from the VLA data, we have used data from the Giant Metrewave Radio Telescope (GMRT) at 350 MHz (Meng et al., in prep). We have recently observed a large region in the Galactic center around SgrB2, reaching an angular resolution of 12.2 ×11.7 (PA = 57 • ) and an rms noise level of 0.05 mJy/beam. In Appendix B we compare the flux of DS at 350 MHz with the VLA data. The spectral index between 350 MHz and 4 GHz (see Fig. B.1) is consistent with the VLA spectral index map shown in Fig. 4, confirming the presence of non-thermal emission in DS.

Radio recombination line emission

In addition to the continuum emission, we also observed RRLs in the whole Sgr B2 region. Due to the low sensitivity of these maps, we have smoothed the images to 8 . At this resolution we are still capable to resolve the structure of the DS region. On the smoothed maps, Gaussian fitting was conducted and peak intensity, centroid velocity, line width and integrated intensity were obtained for all the four stacked RRLs (see Fig. C.1).

The integrated intensity maps of the four RRLs show distinct spatial distribution. At 4.4 GHz, the emission appears relatively diffuse, while at 6.8 GHz, different clumpy structures are visible. At 8.9 GHz the northern part of DS is brighter, while at 10.5 GHz, the southern part of DS has more emission. Such a variation is inconsistent with a simple thermal RRL scenario but suggests the presence of other excitation mechanisms, which will be discussed in Sect. 4.2. The velocity maps of the four RRLs are in agreement for the four images. The eastern lobe of the bubble shows a velocity gradient ranging from velocities about 70 km s -1 in the center down to 55 km s -1 in the outer edge. As shown in the line width maps, the RRLs are typically broad, with values above 30 km s -1 , significantly exceeding the thermal broadening of RRLs in Hii regions with electron temperatures T e = 10 4 K (∆v ≈ 20 km s -1 ). Notably, the line width at the eastern edge of the bubble reaches more than 40 km s -1 . The velocity difference between the eastern edge of DS and the center, together with the increase of line width at the eastern edge, suggests a possible interaction between the expanding Hii region and its surrounding material.

Thermal and non-thermal components in Sgr B2(DS)

The power-law fitting of the continuum emission in Sgr B2(DS) results in a spectral index α that varies from -1.2 to -0.1. The observed range of values of α implies that the centimeter con-tinuum emission in DS is a mixture of thermal and non-thermal contributions. In this section, we dissect the thermal and nonthermal components in the continuum emission and also analyze the properties of the RRL emission.

Disentangling the thermal and non-thermal components

In the following we present two different approaches to disentangle the contributions of the thermal and non-thermal emission in order to better characterize their origin in Sect. 5. The two methods are extrapolating high-frequency emission and fitting the SED with fixed spectral indices.

Extrapolating high-frequency emission

Thermal emission at radio wavelengths is characterized by a relation in which the intensity increases with or is independent of the frequency. On the contrary, the non-thermal emission is characterized by the intensity decreasing with frequency. This suggests that the emission at higher frequencies (corresponding to 11.2 GHz in our dataset) is likely to be dominated by the thermal component, while the emission at lower frequencies (corresponding to 4 GHz) is dominated by the non-thermal component.

In our first approach, we assume that the emission at the highest frequency in our data is dominated by pure thermal (free-free) emission.

The total flux of DS at 11.2 GHz is 0.5 Jy within a diameter of ∼36 , corresponding to a brightness temperature of 14 K. For a typical Hii region temperature of 5 × 10 3 -10 4 K, the optical depth τ ranges from 1.4 × 10 -3 to 2.8 × 10 -3 . Therefore, the free- free emission of DS is optically thin. We use the typical spectral index α = -0.1 of optically thin free-free emission to extrapolate the 11.2 GHz flux density to 4 GHz. The extrapolated thermal component is subtracted from the observed flux density at 4 GHz to get a pure non-thermal component. In Fig. 6 we show the derived thermal and non-thermal components at 4 GHz. The non-thermal emission appears more widespread, while the thermal component appears concentrated in different clumps located along the edge of the bubble. As expected, the Hii region AA has strong thermal emission, while at 4 GHz, the contribution from the non-thermal component drops below the RMS level.

Only in the southeastern region, connecting source AA with the bubble-like structure of DS, we find some presence of a possible contribution of non-thermal emission. The total observed flux at 4 GHz (within the circle highlighted in Fig. 6) is ∼1.5 Jy, for which we determine that ∼60% has a non-thermal origin. We note that since we have assumed that the emission at 11.2 GHz is purely thermal, the contribution of thermal emission at 4 GHz is most likely overestimated and therefore, the non-thermal contribution may be underestimated.

Fitting the SED with fixed spectral indices

In a second approach, we aim at simultaneously determining the contribution of thermal and non-thermal emission. For this, we fit the SED covering the whole range from 4 to 12 GHz with two power-law functions describing each component. For the thermal component, the power-law is S th (ν) ∝ ν -0.1 , assuming that the emission is optically thin free-free. For the non-thermal component, we use S nt (ν) ∝ ν -0.7 (see Hollis et al. 2003;Protheroe et al. 2008;Jones et al. 2011). We fit the observed SED with a linear superposition of these two power-law functions and get the contribution of each component at different frequencies.

As shown in Fig. 7, at 4 GHz, the spatial distribution of the non-thermal component is more widespread compared to the thermal component which is mainly concentrated in the central part of the eastern lobe. This is consistent with the spatial distribution of the two components shown in Sect. 4.1.1. The Hii region AA appears to have pure thermal emission in the center, with a possible contribution of non-thermal emission in the outskirts. The presence of non-thermal emission in the outskirts of Hii regions have been found in a handful of objects (e.g., Garay et al. 1996;Mücke et al. 2002;Veena et al. 2016Veena et al. , 2019)).

From the simultaneous fit of a thermal and non-thermal components, we determine that 90% of the total emission at 4 GHz has a non-thermal origin in Sgr B2(DS). The difference of the relative contribution of thermal and non-thermal components derived from the two methods (see Sect. 4.1.1) can be due to either the assumption of the thermal dominance at 11.2 GHz considered in the model presented in the previous section, or the fixed values of α for the thermal and non-thermal components used in this method. Overall, the results of both methods are in agreement and confirm the presence of extended non-thermal emission in Sgr B2(DS). Further observations at higher and lower frequencies may help to better constrain the properties and distribution of the thermal and non-thermal components in this region.

Stimulated RRLs

As discussed in Sect. 3.3, the RRLs integrated intensity distribution varies significantly among the four stacked frequencies. From Fig. C.1, we can see that the peak intensity of RRLs at the center of the eastern lobe increases monotonically with frequency, which suggests that the RRLs at this part of DS are excited under local thermodynamic equilibrium (LTE) conditions, since under LTE the peak intensity of RRLs obey S ν ∝ ν. However, at the edge of DS, as shown in Fig. 8, the peak intensities of the RRLs exhibit anti-correlation with frequency at 6.8, 8.9, and 10.5 GHz, which suggests that the RRLs in these regions are probably under non-LTE conditions.

We also find that the RRLs are likely excited under non-LTE conditions when comparing their emission with the continuum brightness. Under LTE, the integrated intensity of RRLs relates to optically thin free-free emission as

S L S C ∆v km s -1 = 6.985 × 10 3 ν GHz 1.1 T e K -1.15 , (1)
where S L and S C are the peak flux of RRL and flux of freefree continuum, respectively, and T e is the electron temperature. Adopting T e = 8000 K (see Mehringer et al. 1993), we derive the corresponding free-free continuum level of the four stacked RRLs (RRL-derived free-free emission, hereafter RFE).

The comparison between the RFE and the continuum emission at the four frequencies is shown in Fig. 9. The observed continuum emission is supposed to be a mixture of free-free and synchrotron emission, which should have a higher intensity than the RFE maps. However, the RFE exceeds the observed continuum by a factor of ∼2 at the low frequency end. One possible origin of this excess is that the RRLs are not under LTE conditions but stimulated.

To characterize the stimulation of RRLs, we introduce the RRL peak ratio η. According to theoretical models, the peak intensity of stimulated RRLs and frequency are anti-correlated at high frequencies (Shaver 1978). Therefore, we consider the

where S 10.5 GHz and S 8.9 GHz are the RRL peak intensities, respectively, while the normalization factor 0.85 = 8.9/10.5 result in η = 1 when LTE conditions hold. When the RRLs are stimulated, η < 1. In Fig. 10 we show the η map of DS. From the plot, we see that in the central and western part of DS, η ≈ 1, which indicates that the emission is under LTE. At the eastern edge of DS, we find η < 1, which suggests the presence of stimulated emission. Notably, source AA displays η = 1, meaning that its emission is under LTE, which is consistent with the scenario that source AA is an Hii region with only thermal emission.

We quantify the correlation between non-thermal emission and the stimulation of RRLs, by plotting η against α pixel by pixel (see Fig. 11). A total of 63% of the pixels in DS have α < -0.1, meaning associated with non-thermal emission. Among these non-thermal pixels, 86% (53% of all the pixels in DS) have η < 1 , i.e., the RRLs are stimulated under non-LTE conditions. Out of all the pixels, 58% have η < 1. The vast majority (92%) are associated with non-thermal emission. Therefore, most of the area in DS associated with non-thermal emission is also associated with stimulated RRLs, and similarly, almost all the parts associated with stimulated RRLs show non-thermal emission. The correlation between α and η seen in Fig. 11 is consistent with the scenario that RRLs can be stimulated by non-thermal emission (Shaver 1978).

Origin of thermal and non-thermal emission

We have found that the DS region, located within the envelope of Sgr B2, shows a mix of thermal free-free and non-thermal synchrotron continuum emission at radio wavelengths. This kind of radio continuum emission may have different origins: Hii regions, supernova remnants (SNRs), planetary nebulae (PNe) or extragalactic radio sources. Although the most probable origin is the Hii region, we have explored the other three scenarios. We searched the XMM-Newton catalog for X-ray sources in the region, and found no source associated with DS in the 0.5-12 keV continuum bands (Ponti et al. 2015). The lack of bright X-ray emission suggests that we can exclude a scenario in which DS is a SNR. In a different scenario, PNe can emit in the radio regime, however, their typical size ranges from 0.03 to 0.1 pc, i.e., significantly smaller than the radius of DS (0.5 pc). Additionally, the age of PNe is expected to be between 10 7 and 10 10 yr (e.g., Bressan et al. 1993), much longer than the age of Sgr B2 (estimated to be about 0.74 Myr, see Kruijssen et al. 2015). Also, we estimate the total mass of the ionized gas in DS to be ∼ 500 M , which is much larger than the mass of a typical low-mass star that can generate a PN. Therefore, we exclude the possibility that DS is a PN. The minimum detectable flux density of our images is ∼0.2 mJy, which can be used to estimate the expectable number of extragalactic radio sources. In the entire area of Sgr B2 (7 ×7 ), the expected number of extragalactic radio sources is ∼1, while in the 1 ×1 area of DS, the number is ∼0.03 (see Condon et al. 1998;Anglada et al. 1998). Additionally, the spatial extension of DS (∼1 ) is significantly larger than the typical size of extragalactic radio sources (see Condon et al. 1998). Therefore, the possibility that DS is an extragalactic radio source is excluded. Thus, the only remaining scenario is DS being an Hii region. In this section, we discuss the properties of the central star and the possible mechanisms that can produce the observed non-thermal emission. Fig. 9: Top panels: Free-free continuum maps of the Sgr B2(DS) region derived from the four staked RRLs (RFE as described in Sect. 4.2, following Eq. 1) at 4.4, 6.8, 8.9, and 10.5 GHz.. Bottom panel: Observed continuum emission at 4.4, 6.8, 8.9 and 10.5 GHz. For comparison, the observed continuum maps are also smoothed to 8 , the beam is shown as the shaded circle in the lower left of each plot.

Ionization by a central star

The flux of Lyman continuum photons, ṄLy , needed to ionized an Hii region can be derived as following (Schmiedeke et al. 2016): ṄLy s -1 = 8.9 × 10 40 S ν Jy

ν GHz 0.1 T e 10 4 K -0.45 D pc 2 , (3)
where S ν is the flux at frequency ν, T e is the electron temperature, and D is the distance to the source. The total flux of the free-free emission in DS is 0.5 Jy at 11.2 GHz, assuming that the continuum emission at this frequency is pure thermal. For a typical Hii region, the electron temperature ranges from 5000 to 10000 K (Zuckerman et al. 1967). Taking D = 8.34 kpc for Sgr B2, we obtain a flux of Lyman continuum photons of 4-6 ×10 48 s -1 . Such a Lyman continuum flux corresponds to an O7 ZAMS star (see Table II of Panagia 1973). We have searched the Spitzer infrared images (e.g., Ramírez et al. 2008) for a possible infrared counterpart of the star ionizing the Hii region in DS, but found no clear candidate. Also the young stellar object catalog of (Yusef-Zadeh et al. 2009) does not show any infrared source in the center of DS. One possible reason for the lack of a detected infrared source in DS is the high extinction. Assuming a column density of molecular hydrogen of 10 24 cm -2 in DS (Schmiedeke et al. 2016;Ginsburg et al. 2018), the extinction at 5 µm is ∼ 20, based on the grain model by Li & Draine (2001).

For an O7 star at a distance of 8.34 kpc, this extinction result in an apparent magnitude ∼ 30, which is non-detectable in current infrared images.

Non-thermal emission origin

The interaction between cosmic-ray nuclei and the interstellar medium produces secondary electrons (and positrons). These secondary electrons can be a source of synchrotron emission like the one detected in Sgr B2(DS). Protheroe et al. (2008) analyze diffusion models and indicate that for secondary electrons to radiate synchrotron emission at GHz frequencies in ∼1 mG magnetic field, cosmic-ray particles should have multi-GeV energies, whereas such cosmic-ray particles cannot penetrate into the cloud due to the high gas density in Sgr B2 (see, e.g., Padovani et al. 2009Padovani et al. , 2013)). Thus, the penetration of cosmic rays from outside (e.g., from the nearby Galactic center) is excluded as a possible origin to the non-thermal emission observed in Sgr B2(DS). However, the presence of synchrotron emission is the fingerprint of relativistic electrons. Padovani et al. (2015Padovani et al. ( , 2016) ) showed that particle acceleration can take place at a shock loca- tion in protostars according to the first-order Fermi acceleration mechanism (also known as diffusive shock acceleration). In fact, local thermal particles are accelerated up to relativistic energies while crossing back and forth a shock surface thanks to magnetic fluctuations around the shock itself.

In a companion paper (Padovani et al. 2019), we present an extension to this theory applied to Hii regions. We demonstrate that electrons can be efficiently accelerated under the hypothesis that the flow velocity in the reference frame of an expanding Hii region hitting denser material is sufficiently high (> 35 km s -1 ) to switch on particle acceleration. The main parameters of the model are the magnetic field strength, the the flow velocity in the shock reference frame (hereafter velocity), the temperature, the volume density, and the ionization fraction. Assuming a complete ionized medium and a temperature of 8000 K (see Sect. 4.2), we studied the parameter space volume density-magnetic field strength and computed the emerging shock-accelerated electron flux for a shock radius of 0.36 pc, which is the average distance between the center of DS and the inner surface of the synchrotron emitting region. This electron flux is used to compute the synchrotron emissivity and then the flux density between 4 and 12 GHz for a beam of 4 and an average emitting region size3 of 0.85 pc. We performed a χ 2 test 3) by the First-order Fermi acceleration model through a χ 2 test using the model described in Padovani et al. (2019). The model also generates the modeled spectral index (α mod ) map which is consistent with the observed α obs map.

to find the values of the velocity (U), density (n), and magnetic field strength (B) that best reproduce the observations. For the flux densities we obtain an average accuracy of about 20% assuming 33 U 50 km s -1 , 10 4 n/cm -3 9 × 10 4 , and 0.3 B/mG 4. From the model, we derive maps of the velocity, magnetic field strength and density in the Sgr B2(DS) region (see Fig. 12). The magnetic field strength is found to be in the range 0.5-2 mG as reported by Crutcher et al. (1996). The model also gives an α map of DS, which is in agreement with the observed α map (Fig. 5). It is interesting to notice that the velocity is lower where the density is higher. One reason might be that the flow is moving faster towards the east-west direction, which is the direction of minimum resistance. In contrast, the velocity is lower toward the northern direction where the density is higher. Although the RRL line width cannot be directly used to determine the shock speed, we note that it is significantly larger than the typical expansion velocity of an Hii region (10 km s -1 Draine 2011), which points to an additional mechanism of acceleration. One possible origin of the high velocity is the presence of a stellar wind driven by the O7 star that likely ionizes the Hii region (see e.g., Veena et al. 2016;Pereira et al. 2016;Kiminki et al. 2017). We note that velocities of about 35 km s -1 are in agreement with what is obtained by simulations of Hii regions of O and B stars driving strong stellar winds (Steggles et al. 2017).

It is important to remark that the flux densities (from which the spectral index is derived) estimated by our model are based on the assumption that all the parameters (magnetic field strength, velocity, temperature, density, and ionization fraction) are constant along the line of sight. For a proper modeling, one should account for the spatial variations of these quantities. However, the fact that model results are already fairly comparable with observations is an indication that shock-accelerated electrons may represent the key to explain the synchrotron emission in Sgr B2(DS). Previous studies (e.g., Yusef-Zadeh et al. 2007a,b, 2013, 2016) indicate that cosmic-ray induced nonthermal emission is common in the Galactic center region. On large scales, where the density is not as high as in the envelope of SgrB2, cosmic rays may play a dominant role in the production of relativistic electrons, while in extremely dense regions (and smaller scales, like Sgr B2(DS)) Hii region shocks may be a dominant source of non-thermal emission.

Summary

We have observed the giant molecular cloud Sgr B2 with the VLA in its D and CnB configurations. We obtained continuum and RRL maps in the frequency range 4-12 GHz, covering a spatial extent of about 20 × 20 , with a beam size of 4 . We have focused our study on the Sgr B2(DS) region, located in the southern area of the Sgr B2 envelope. Our main results are:

-Sgr B2(DS) is bright at radio wavelengths, with intensities in the range 10-50 mJy/beam, within the observed 4 to 12 GHz frequency range. At 4 resolution, DS appears as bubble-like Hii region with a diameter of about 1.5 pc, powered by an O7 star, and surrounded by dense gas and a series of dense cores distributed in an arc structure around it. -We find that the total flux of DS decreases from 1.9 Jy at 4 GHz down to 0.5 Jy at 12 GHz. Spectral analysis shows that the spectral index of DS varies from -0.4 to -1.2, suggesting the presence of non-thermal emission, in addition to the thermal emission of the Hii region. We decompose the thermal and non-thermal components in DS, and find that the thermal emission is clumpy and concentrated in the eastern lobe of the bubble-like structure, while the non-thermal emission appears widespread over the whole region. -The emission of RRLs in DS show a central velocity varying from 70 km s -1 in the center to 55 km s -1 in the outer edge.

The line widths of the RRLs range from 30 to 40 km s -1 .

From the RRLs integrated intensity, we derive the corresponding thermal continuum emission, which exceeds the observed continuum by a factor of about 2 at the low frequency end. In addition, the RRLs intensity does not follow the S ν ∝ ν relation expected for RRLs under LTE conditions, but drops at high frequencies likely due to non-LTE effects. We find a correlation between the presence of nonthermal emission (i.e., negative spectral indices) and those regions where RRLs are excited under non-LTE conditions. This suggests that RRLs in Sgr B2(DS) are possibly stimulated by synchrotron emission. -We modeled the observed synchrotron emission, and found that relativistic electrons can be produced via first-order Fermi acceleration, which is triggered by the interaction between the expanding Hii region and the denser surrounding From the model, we derive maps of the flow velocity in the shock reference frame, magnetic field strength and density in the Sgr B2(DS) region. Velocities are found to be between 35 and 50 km s -1 as found in simulations of cometary Hii regions of O and B stars driving strong stellar winds. The magnetic field strength is found to be in the range 0.3-4 mG as reported by Crutcher et al. (1996) and the density in the range of 1 -9 × 10 4 cm -3 .

In this work we have characterized the non-thermal emission observed towards the Hii region Sgr B2(DS), and we have proposed that locally accelerated particles, under certain conditions of density, magnetic field, and shock velocity, can be the origin of the non-thermal emission. This result raises the question whether other Hii regions may also show non-thermal emission. A handful of Hii regions have been reported to have negative spectral indices suggestive of a non-thermal component (e.g., Garay et al. 1996;Veena et al. 2016Veena et al. , 2019)), indicating that the presence of non-thermal emission in Hii regions may be a more common event than previously thought. The advent of extremely sensitive facilities at radio wavelengths such as SKA (Square Kilometer Array) or the ngVLA (next generation Very Large Array) may trigger the systematic search for Hii regions with nonthermal emission, thus allowing a study of the conditions under which non-thermal emission exists.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Continuum images of Sgr B2 in C (panel a) and X (panel b) bands. Relevant regions are marked with their names (see Mehringer et al. 1993). The dashed boxes mark the region of DS. The synthesized beam of the C and X images are 2.7 ×2.5 and 1.82 ×1.53 , respectively

Fig. 3 :

 3 Fig. 3: Channel maps of the Sgr B2(DS) region. All the 12 maps have been produced considering the same uv range limited to 0.6-50 kλ, and have been convolved to a circular beam of 4 . The synthesized beam is shown as a yellow circle in the bottom left panel.

Fig. 4 :Fig. 5 :

 45 Fig. 4: The spectral index (α, panel a) and its uncertainty, (panel b) throughout the whole Sgr B2 region. The regions marked in both panels correspond to those regions also labeled in Fig. 1. See Sect. 3.2 for details in the calculation of the spectral index.

FFig. 6 :Fig. 7 :

 67 Fig. 6: Spatial distribution of synchrotron (or non-thermal, left panel), mixed (central panel, corresponding to observed image) and free-free (or thermal, right panel) components of DS at 4 GHz as derived from method described in Sect. 4.1.1. The thermal emission map is extrapolated from 12 GHz assuming that at this frequency the emission is originated by a thermal component with a spectral index α = -0.1. The non-thermal emission is obtained by subtracting the extrapolated thermal emission from the mixed emission. The contours are the same as in Fig. 5. All panels have a circular beam of 4 (shown in the bottom left corner of each panel). The dashed circles indicate the region in which we calculate the flux density of DS (see Sect. 4.1).

Fig. 8 :

 8 Fig. 8: Right: Spectral index map of the Sgr B2(DS) as shown in the central panel of Fig. 5. The map of the uncertainty of spectral index is as shown in the right panel of Fig. 4. Left: Peak intensities of the four stacked RRLs (see Table1) at four selected positions. Position (a) corresponds to the well-known Hii region AA, while the other positions have been selected to probe regions with negative spectral indices. The dashed line in each panels marks the scenario for LTE under which the flux is proportional to the frequency.

FFig. 10 :

 10 Fig. 10: Map of η (see Sect. 4.2) in DS. The X band continuum emission is overlaid as contours. The synthesized beam, corresponding to 8 , is plotted as a dark circle at the lower left corner. The dashed circle indicates the area in which pixels are taken into account for Fig. 11.

Fig. 11 :

 11 Fig.11: Plot of η against spectral index (α). Only the pixels that are in the area marked by the dashed circle in Fig.10are taken into account, thus the influence of source AA is excluded. The criterion of non-thermal emission, α = -0.1, is shown with a vertical dashed line, while the criterion of stimulated RRL emission, η = 1, is shown as a horizontal dashed line. These two orthogonal lines divide the whole α-η space into four quadrants. For each quadrant, the number of pixels within it as a percentage of the number of all the pixels is written. The points with error bars crossing the α = -0.1 or η = 1 lines are colored gray and neglected in the statistics.

Fig. 12 :

 12 Fig.12: Maps of shock velocity (U), volume density (n), and magnetic field strength (B) of DS that reproduce the observed flux density maps at the 12 frequencies (see Fig.3) by the First-order Fermi acceleration model through a χ 2 test using the model described inPadovani et al. (2019). The model also generates the modeled spectral index (α mod ) map which is consistent with the observed α obs map.

F.

  Meng et al.: Non-thermal emission in the envelope of Sgr B2 material. The model, which is presented in a companion paper ofPadovani et al. (2019), reproduces the observed flux density and spectral index.

  Fig. A.1: Top panel: Continuum emission image of Sgr B2 at 4 GHz. The dashed line indicates the trajectory along which we get the intensity profile. Bottom panels: Intensity profiles of the twelve channel maps. The uv limit of the twelve images has been set to 0.6-50 kλ and the images have been convolved to a circular beam of 4 . The dashed line indicates the zero intensity level.

FigFig

  Fig. B.1: Right panel: Spectral index image obtained from the 350 MHz GMRT emission map (Meng et al., in prep.) and the 4 GHz VLA continuum image. Both images have been obtained after limiting the short uv to 0.6 kλ, and convolving the images to a common beam of 12.2 ×11.7 (PA = 57 • ), which corresponds to ∼ 25 kλ. Left panels: Two examples of SEDs for two selected positions within DS. The observed fluxes are shown with green symbols, all of them measured from images with same spatial filtering and synthesized beam. For illustration, the thermal (blue dashed curve), non-thermal (red dashed line) and the mixed (black solid curve) components are marked. The curves are just qualitative descriptions of the possible contribution of the thermal and non-thermal components, and do not aim at fitting the data.

Table 1 :

 1 Observed and stacked RRLs.

	RRL	ν 0 (MHz)	Stacked (ν) a
	H115α	4268.14
	H114α H113α	4380.95 4497.78	RRL 4.4 GHz
	H112α	4618.79
	H100α	6478.76
	H99α H98α	6676.08 6881.49	RRL 6.8 GHz
	H97α	7095.41
	H92α	8309.38
	H91α	8584.82
	H90α	8872.57	RRL 8.9 GHz
	H89α	9173.32
	H88α	9487.82
	H87α	9816.86
	H86α	10161.30
	H85α	10522.04	RRL 10.5 GHz
	H84α	10900.06
	H83α	11296.41

A new distance to the Galactic center has been measured to be 8.127 ± 0.031 kpc (Gravity Collaboration et al.

2018). For consistency with the paper published within the same series of studies of Sgr B2, we use the distance reported byReid et al. (2014).

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Article number, page 4 of 15 F. Meng et al.: Non-thermal emission in the envelope of Sgr B2

The average size, L, of a shell-shaped emitting region with inner and outer radii, R in and R out , respectively, is L = (π/2)R out (1 -R 2 in /R 2 out ). For DS we assume R in = 0.36 pc and R out = 0.72 pc.

Acknowledgements. We thank the referee for his/her comments during the reviewing process.. FM, ASM, PS, ASchw research is carried out within the Collaborative Research Centre 956, sub-projects A6 and C3, funded by the Deutsche Forschungsgemeinschaft (DFG) -project ID 184018867. MP acknowledges funding from the European Unions Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 664931. This research made use of Astropy, 4 a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018).

Appendix A: Intensity profiles

In Fig . A.1 we plot the intensity profiles of the twelve channel maps across the Sgr B2(DS) region. The profiles are taken along the RA direction from (17:47:05, -28:25:10) to (17:47:38, -28:25:10). Since we set uv limit to 0.6-50 kλ, the large scale structures are filtered out, resulting in the background intensity to be ∼ 0.

Appendix B: Combining GMRT and VLA

In this section we make use of recently observed GMRT observations to confirm the negative spectral indices measured towards the Sgr B2(DS) region. The GMRT observations were taken under the project number (31_021) and will be presented in detail in a forthcoming publication. In summary, we used the GMRT to observe the Galactic center at 350 MHz, with a total bandwidth of 32 MHz. The observations were centered at the position of Sgr B2. For comparison with the VLA data, we have selected data in the uv range 0.6-25 kλ. The short spacing is set to 0.6 kλ to match the same scales recovered in the VLA data. The synthesized beam of the GMRT image is 12 . 2 × 11 . 7 (PA = 57 • ). We have convolved the VLA images to the same resolution of the GMRT data, which results in the upper limit of the uv range for both the GMRT and VLA images to be ∼25 kλ.

We obtain the spectral index in DS by interpolating the flux density between the 350 MHz and 12 GHz images. As shown in the right panel of Fig. B.1, the spectral index in DS is negative and in agreement with the spectral index values determined using the VLA data (cf. Fig. 4). In the left panels of Fig. B.1 we show the SED of the 13 channel images for two positions in DS. The SEDs in these positions can be describe by a combination of thermal and non-thermal emission, in agreement with the analysis based on the VLA data. However, the relatively coarse angular resolution of the GMRT images makes the study of the small-scale structure of the DS bubble not possible. Therefore, in this paper, we focus on the analysis of the VLA images and the sub-pc structures in DS. A complete analysis of the GMRT data will be used to study the large-scale structure of the Galactic center region in a forthcoming paper (Meng et al., in prep.).

Appendix C: The RRL maps

In this section we show the maps of fit parameters of the four stacked RRLs (see Table 1). The fit parameters include integrated intensity, centroid velocity, line width, and peak intensity (see Fig.