THz Transmission meets Untrusted UAV-Relaying; Trajectory and Communication Co-design for Secrecy Energy Efﬁciency Maximization

—Unmanned aerial vehicles (UAVs) and Terahertz (THz) technology are envisioned to play paramount roles in next-generation wireless communications. Hence, this paper presents a novel secure UAV-assisted mobile relaying system operating at THz bands for data acquisition from multiple ground user equip- ments towards a destination. We assume that the UAV-mounted relay may act, besides providing relaying services, as a potential adversary called the untrusted UAV relay. To safeguard end-to-end communications, we present a secure two-phase transmission strategy with cooperative jamming. Then, we formulate an optimization problem in terms of a new measure − secrecy energy efﬁciency (SEE), deﬁned as the ratio of achievable average secrecy rate to average system power consumption, which enables us to obtain the best possible security level while taking UAV’s inherent ﬂight power limitation into account. This optimization problem leads to a joint design of key system parameters, including UAV’s trajectory and velocity, communication scheduling, and power allocations. Since the formulated problem is a mixed-integer non-convex optimization and computationally intractable, we propose alternative algorithms to solve it efﬁciently via greedy/sequential block coordinated descent, successive convex approximation, and non-linear fractional programming techniques. Numerical results demonstrate signiﬁcant SEE performance improvement of our designs when compared to other known benchmarks.


I INTRODUCTION
T HE unmanned aerial vehicle (UAV) has recently been recognised as one of the major technological breakthroughs to be pervasively applied in 5G-and-beyond wireless communication networks supporting massive machine-type communications, internet of things (IoT), and artificial intelligent (AI) empowered communications [1]- [3]. Thanks to the unique characteristics of agility, on-demand swift deployment, versatility, and channel superiority amongst the other potentialities, UAV-aided wireless communications have recently attracted a great deal of research [4]- [7]. Despite numerous advantages, the open nature of air-ground (AG) links inevitably makes such systems vulnerable to malicious attacks such as eavesdropping. Accordingly, the security and confidentiality of such promising wireless communication systems are of utmost concern and undeniable requirements. To protect the confidentiality of UAV communications against hostile entities, one promising technique is the physical layer security (PLS) that uses the characteristics of wireless channels and applies communication techniques to combat attacks without complex encryption. A number of works have found leveraging the PLS in UAV-aided communications plausibly effective [8]- [22]. For example, PLS has been exploited in a wireless-powered UAV-relay system to combat eavesdropping via maximizing secrecy rate by a joint design of UAV's position and resource allocation [10]. Other efforts were made to maximize the average secrecy rate (ASR) via joint trajectory and communication design for UAV-standalone wireless system [11]- [13], for double-UAV with external jamming [14]- [16], and for secure UAV-relaying scenarios [17]- [22]. The majority of previous research has deemed the UAV to be a fully authorized and legitimate communication node in UAV-assisted relaying applications. However, when the UAV behaves as an untrusted relay, which is called untrusted UAV-relay (UUR), with the capability of information eavesdropping while assisting endto-end communications (see [23], [24]), the system design becomes quite challenging and entirely different from the existing body of research.
Further, energy efficiency is another imperative need for UAV-aided communications due to UAVs' inherent constraints on size, weight, and power (SWAP). Typically, the small-scale rotary-wing UAVs are powered via limited on-bored batteries, leading to a restrictive operational lifetime, which undoubtedly impacts their overall system performance. Nonetheless, UAVs' flight endurance, if properly designed, can be enhanced to a considerable extent [25]. Several works have studied the secrecy performance of UAV-aided systems considering the propulsion energy consumption constraint [26]- [29]. In [26], the authors have investigated ASR maximization for a cooperative dual-UAV secure data collection with propulsion energy limitation. Exploring the problem of secrecy energy efficiency (SEE) maximization for UAV-aided wireless systems is another research path [27]- [29]. The authors have designed both trajectory and resource allocation for the energy-efficient secure UAV communication system with the help of a multiantenna UAV-jammer in [27]. Some appropriate system designs have been conducted for the SEE improvement of a single UAV-relay system [28], and a UAV-swarm multi-hop relaying scenario [29]. It is worth pointing out that all the aforementioned designs have only aimed to combat external terrestrial eavesdroppers.
On the other hand, owing to the ultra-broad bandwidth at the terahertz (THz) frequency range (0.1-10 THz), THz transmission has been acknowledged as a promising technology capable of catering an explosive growth of user demand of higher mobile traffic for future wireless systems. However, THz links incur severe path loss and high susceptibility to environmental blockage, and molecular absorption [1], [5], [30], which limit signal propagation distance and coverage range. To overcome the hindrances, one possible solution might be to explore UAVaided communications in THz links. Notably, in the context of THz-UAV systems, few initial research studies have thus far been conducted. The coverage probability of the UAV-THz downlink communications was analyzed in [31], while [32] has explored a similar non-security scenario with a focus on minimizing communication delay by a joint design of the UAV's location and power control. When it comes to security issues of such high-frequency systems, despite the widelyassumed improved resiliency against eavesdropping of THz links, the authors of [33] have characterized the possibility of eavesdropping attacks for such systems. Needless to mention that even considering negligible information leakage towards the external malicious eavesdroppers through THz transmissions, the scenarios involving untrusted relays, particularly the UUR systems, may still be vulnerable to eavesdropping. The appropriate design for such systems has yet to be understood; therefore, one needs to design novel frameworks to enable the efficient deployment of THz-UUR wireless systems.

I-A Our contributions
To the best of our knowledge, this is the first work addressing the energy-efficient secure design of a THz-UUR wireless communication system to guarantee confidentiality and transmission secrecy with the least system power consumption. Our detailed contributions are summarized below.
• We present an UUR-enabled wireless communication system for data collection from multiple ground user equipments (UEs) towards the base station (BS) over THz-based AG links. We adopt a secure two-phase transmission strategy using destination-assisted cooperative jamming to improve security. • Then, we formulate a maximin optimization problem in terms of a new measure secrecy energy efficiency (SEE), defined as the ratio of achievable ASR to average system power consumption. This optimization problem leads to a joint design of key system parameters, including UUR's trajectory and velocity, communication scheduling, and network transmission power allocations. • Since the optimization problem is originally intractable due to non-convexity, we decompose it into four subproblems and then solve each via successive convex approximation (SCA) or Dinkelbach fractional programming techniques. Further, we propose two computationally efficient algorithms according to the sequential and maximum improvement (MI) based block coordinate descent (BCD) approaches with guaranteed convergence to at least a suboptimal solution. We also conduct computational and complexity analysis and show that our solution can be obtained in polynomial time order, making it applicable to the energy-hungry UAV-based scenarios. • We conduct extensive simulations to verify the analyses and demonstrate the effectiveness of our proposed designs in terms of SEE compared to some other benchmarks, without communication resource allocation design or trajectory and velocity optimization and ignoring flight power consumption. We also investigate the impact of some fundamental setting parameters such as the flight mission time and the molecular adsorption factor on the overall system secrecy performance.
The rest of the paper is organized as follows. Section II introduces system model and formulates the problem of interest. In Section III, we present efficient iterative algorithms to solve the optimization problem, followed by numerical results and discussions given in Section IV. Finally, we draw conclusions in Section V.

II SYSTEM MODEL AND PROBLEM FORMULATION
We consider a UAV-enabled wireless communication system for data collection from a set of K ground UEs towards a BS via a UAV-assisted mobile amplify-and-forward (AF) relay, as shown in Fig. 1. Here we assume that there are no reliable direct links from UEs to BS (see [24], [34] and references therein), and all nodes are equipped with a single antenna, operating in half-duplex mode. Therefore, a UAVrelay is employed to assist end-to-end communications [18]; nonetheless, the UAV-relay may not be fully authorized to access collected confidential information and may conduct malicious eavesdropping, i.e., an UUR [23]. Thus secure data transmission is in demand.
Without loss of generality, we consider a 3D Cartesian coordinate system, where the BS's horizontal coordinate is located at the origin q b = [0, 0] ∈ R 1×2 , and the ground UEs with horizontal coordinates q k = [x k , y k ] ∈ R 1×2 for ∀k ∈ K, where K = {1, 2, · · · , K}, are randomly distributed in a circular annulus region with the inner radius R 1 and outer radius R 2 and the coordinates are assumed to be known in prior. Here, R 1 is considered to be the largest distance at which a reliable uplink transmission can be obtained, while beyond R 1 in our case implies no direct link between US and BS. Further, R 2 indicates the boundary of the permitted flying region for the UAV to provide communication service.
We also consider that UAV flies from and back to the same specific point over the region of interest for a duration of T seconds in order to provide relaying services to all UEs with fairness. This specific point may refer to the checkup point wherein the UAV gets recharged and physically examined to maintain its service. Assuming that UAV flies at a fixed altitude H meters whose instantaneous horizontal coordinate and velocity is represented by and v(t) dt , respectively, where 0 ≤ t ≤ T . For the ease of analysis, we adopt the time-slotted system such that the flight duration T is equally discretized into N sufficiently small time slots of duration δ t

II-A Channel model
We assume that the AG links are over THz channels which also follow reciprocity. Therefore, assuming that, at each time slot n, the channel state information is regarded static due to adequately small δ t , we adopt the line-of-sight (LoS) dominant THz channel power gain model [35] between the UUR and any UE k ∈ K as where d ku denotes the Euclidean distance between the UUR and the k-th UE, given by Note that the multiplicative term exp(−a f d ku ) in (1) is the indication of excessive path loss of THz links due to water vapor molecular absorption effect, wherein a f ∆ = a(f ) is the frequency dependent adsorption factor, β 0 ∆ = ( C 4πf ) 2 denotes the reference channel power gain at unit distance, wherein C is the speed of light, f is the operation frequency. Likewise, the THz channel power gain between the UUR and the BS can be written as

II-B Constraints on user scheduling, power, UAV's mobility
We adopt the time division multiple access (TDMA) protocol for multiuser relaying services, wherein UUR serves at most one enabled UE at n-th time slot, while the other ground UEs keep silence. Therefore, letting ζ k [n] be a binary user scheduling variable for UE k ∈ K at time slot n ∈ N , we have the user scheduling constraints C2 : where ζ k [n] = 1 if user k is scheduled at time slot n, and ζ k [n] = 0, otherwise. Further, the transmit powers of the UUR, the BS, and k-th user in time slot n, denoted respectively as , and p k [n], are generally subject to average and peak transmit powers given as where sets {p ave u , p ave b , p ave k , ∀k} and {p max u , p max b , p max k , ∀k} represent the corresponding average and maximum power constraints of the network nodes.
The mechanical power consumption of the energy-limited UAV due to high demand of propulsion energy for aerial operation can be approximately given by [25] wherein v[n] is the UAV's instantaneous velocity at time slot n, P 0 and P i are two constants representing UAV's blade profile power and induced power in hovering mode, respectively, Ω u and R u are the UAV's blade angular velocity in Radian per second (rad/s) and its rotor radius in meter (m), d 0 , ρ, s, and A indicate the unit-less fuselage drag ratio, air density in kg/m 3 , rotor solidity, and rotor disk area in m 2 , respectively. Further, the average rotor induced velocity in hovering is shown as v 0 . And we have the average flight power consumption constraint as whereinP lim indicates the UAV's average propulsion power budget, which is proportional to the UAV's onboard battery capacity. Therefore, it should be required that the total consumed propulsion energy by the UAV over N time slots be less than such limit in order for network functioning. Further, the considered scenario should be subject to UAV's mobility constraints in terms of initial and final locations for cyclic path, in-flight maximum displacement per time slot for satisfying channel invariant assumption, and permitted flying region as wherein q I indicates UAV's initial and final location per flight, v max u and a max u are the UAV's maximum speed and acceleration, respectively.

II-C Secure transmission strategy, problem formulation
For the purpose of wireless security, we adopt a secure twophase transmission strategy with destination-assisted cooperative jamming (DBCJ) technique similar to [11], [18], [24]. In the first phase, at each time slot n, the scheduled UE sends confidential information to UUR, and simultaneously the BS jams UUR. In the second phase, UUR forwards the received signals using AF relaying protocol to the BS.
Under such setting, given the equally shared communication bandwidth B Hz, the achievable end-to-end instantaneous data rate in bits-per-second (bps) from the k-th UE at time slot n is Then the UUR may overhear the confidential information with an achievable wiretap secrecy rate per Hz at time slot n as indicate the power of additive white Gaussian noise (AWGN) at the receivers, wherein σ 2 u and σ 2 b , which are assumed equal for simplicity, denote the power spectral density (PSD) at the UUR and the BS.
We adopt the ASR as one of the key secrecy metrics and the ASR of k-th UE at time slot N is wherein (x) + ∆ = max{x, 0}, and the ratio 1 2 is due to the fact that secure transmission is done in two phases of equal duration at each time slot. The achievable average information bits can securely be exchanged between k-th UE and BS is . To fully exploit the capability of aerial platforms for communication, the limited energy resource must be considered in system design. In practice, the UAV's propulsion power consumption is much higher than those used for UEs' signal transmission, BS's jamming, and receiver processing. Hence, we approximate the network's total power consumption mainly from UAV's propulsion only. Consequently, for the secrecy metric, we define secrecy energy efficiency (SEE) of the proposed scheme for the k-th UE as the ratio of the achievable ASR to the approximated average system power consumption, , bits/Joule (17) wherein the user scheduling set ζ ζ ζ = {ζ k [n], ∀n, k}, UAV's location and velocity set Q = {q[n], v[n], ∀n}, and network transmit power set , ∀n}} are the involving parameters.
Remark 1. It is worth pointing out, for later analysis, that we use normalized metrics, i.e., the numerator and denominator of (17) divided by B andP lim , respectively, to balance well numerical values and both metrics in SEE.
To design the network to obtain the best performance among UEs and provide fairness support to all UEs given UAV's stringent on-board battery, we maximize UEs' minimum SEE by We note that the problem (P) is a mixed-integer non-convex optimization problem, which is too hard to solve optimally. The non-convexity is mainly due to the non-concave objective function with respect to (w.r.t) the optimization variables, and also having the non-smoothness operator (·) + and the nonconvex constraints (C1), (C3), and (C9). To make it tractable, we first remove the operator (·) + from the numerator of the objective function, since the value of the objective function should be non-negative at the optimal point; otherwise, one can set, e.g., P k = 0, ∀k and get zero SEE performance without modifying the original problem. Nonetheless, having at least a differentiable objective function, the problem is still nonconvex, thereby no standard approach to solve it efficiently. To remedy this issue, we propose some computationally efficient algorithms to iteratively solve a sequence of approximated convex sub-problems by adopting several techniques such as block coordinated descent (BCD), successive convex approximation (SCA), and nonlinear fractional Dinkelbach programming, discussed below.

III PROPOSED ITERATIVE SOLUTION
In this section, we split the problem (P) into four subproblems with different blocks of variables, then solve each block by block, while keeping the other blocks unchanged. Specifically, we delve into solving the joint user scheduling and transmit power optimization sub-problem to optimize (ζ ζ ζ, P k ), relaying and jamming power optimization subproblems to improve P u and P b , and lastly, the joint trajectory and velocity optimization subproblem to optimize Q. Then, the overall algorithms to iteratively attain the approximate solution of (18) will be given.

III-A Joint user scheduling and transmit power optimization
First we relax binary variables ζ ζ ζ into continuous real-valued setζ ζ ζ = {ζ k [n], ∀k, n}. We define the auxiliary variables , ∀k, n}. Now, introducing a slack variable ψ and given the local point in the l-th iteration (P , the corresponding relaxed sub-problem can be represented as where Note that the constraint (19a) should be met with equality at the optimal point; otherwise, the value of the objective function in problem (P1) can still be increased by increasing ψ, which violates the optimality. The sub-problem (P1) is still nonconvex due to non-convexity of the constraint (19a) and for general N , it is indeed NP-hard. Therefore, we cannot solve it efficiently. To handle (19a), we first present Lemma 1 below.
Lemma 1. Let's define the bivariant functions Z 1 (x, y; a, b) x ln(1+ ay y+bx ) and Z 2 (x, y; c) ∆ = x ln(1+ cy x ) over the domain x, y > 0 with the positive constants, i.e., a, b, c > 0. Both Z 1 and Z 2 are jointly concave w.r.t the variables x and y. Additionally, the inequality below near the given point (x 0 , y 0 ) always holds with tightness: Proof. Please see Appendix A.
Using Lemma 1, it can be identified that both Terms I and II in (19a) are concave w.r.t the optimization variablesζ ζ ζ and P k , since the summation operator preserves the convexity. The non-convexity of the left-hand-side (LHS) expression is in the form of concave-minus-concave. Then using (20) and applying the SCA technique, we approximate the non-convex constraint with the corresponding approximate convex one at each iteration. Given the local point (P ) in the l-th iteration, we can reformulate (P1) as follows.
Since the reformulated problem (P1.1) is convex w.r.t the optimization variables {ψ,ζ ζ ζ,P k }, it can be solved efficiently via CVX using the interior-point method [36]. Having solved this subproblem, we can then obtain the optimized value of , ∀k, n}. Further, once the solution of overall algorithm is obtained, we can reconstruct the corresponding binary solution of ζ ζ ζ, according to the method in [37], or using ζ ζ The formulated convex optimization model given in (P1.1), though being convex, cannot be directly accepted by CVX, as it does not follow the disciplined convex programming (DCP) ruleset required. Given that the relative entropy function E rel (x, y) = x log( x y ), x, y > 0 is convex and accepted by CVX, we can rewrite concave function Z 1 (x, y; a, b) (or the equivalent expression in the constraint (19a)), as where the equality (a) follows from the following relations between different form of logarithmic functions and the convex relative entropy function given by x ln 1 + wherein (23) and (24) are jointly concave and convex w.r.t the joint variables (x, y) over x, y > 0, respectively.

III-B Relaying power optimization
The corresponding sub-problem for optimizing UUR's relaying power can be rewritten, introducing the slack variable ψ, as and ∀k, n Note that sub-problem (P2) is a convex optimization problem due to having an affine objective function and all convex constraints, following from Lemma 2 introduced below.
f 2 (x) is concave subject to the condition ad ≥ bc, following from the fact that the function ln(1 + qx), q ≥ 0, x > 0 is concave w.r.t x, whose extended-value extension is nondecreasing and h(x) = − 1 x is also concave; therefore, (f • g)(x) is concave. Note that the last equality of (26) represents the understandable reformulation of the function f 2 (x; a, b, c, d) by the CVX optimization toolbox. We also stress that for any given point x 0 , there is a unique convex function f lb 2 (x; x 0 , a, b, c, d) defined as . (27) such that f lb 2 (x; x 0 , a, b, c, d) serves as a global lower-bound of f 2 (x), i.e., f 2 (x) ≥ f lb 2 (x; x 0 , a, b, c, d) [38].
Consequently, one can solve subproblem (P2) efficiently using CVX. Here, we have (N + 1) optimization variables and (N + K + 1) convex constraints. Assuming the convergence accuracy of interior-point algorithm employed for solving this convex problem with logarithmic cone is ε 2 , the complexity cost of solving sub-problem (P2) can be obtained .

III-C Jamming power optimization
Keeping the other variables unchanged and taking the slack variable ψ, the BS's jamming power optimization sub-problem is given as . Notice that subproblem (P3) is non-convex due to non-convex constraint (28b), which is in the form of convex-minus-convex according to [15,Lemma 1]. Therefore, we apply SCA such that for a given local point P in l-th iteration. We approximate the first convex term with the global underestimator concave expression and obtain the convex reformulation as Since subproblem (P3.1) is convex, we can solve it efficiently using CVX. Here, we have N + 1 optimization variables and (N + K + 1) convex constraints. Assuming the accuracy of SCA algorithm for solving this problem is ε 3 , the complexity of solving approximated sub-problem (P3.1) can, therefore, be represented as O (N + 1) 2 (N + K + 1) 1.5 log 2 ( 1 ε3 ) .

III-D Joint trajectory and velocity optimization
Now, we optimize the trajectory q and velocity v of the UUR while keeping the transmit power allocation and user scheduling sets (P, ζ ζ ζ) fixed. Therefore, the corresponding subproblem can be given as . In order to solve subproblem (P4), we should maximize every single fractional terms of , ∀k subject to the given constraint (31b). In light of this, let λ be the maximum SEE of sub-problem (P4) with solution set (q , v ) given by wherein F represents the feasible set spanned by the constraint (31b). Applying nonlinear fractional Dinkelbach programming theory [39], the objective function of problem (P4) can be equivalently transformed into a subtractive version such that the optimal value of λ can be achieved iff Thus, we can optimize the equivalent problem to obtain the optimal solution of Q, via solving the reformulated problem as showing the value of λ in the m-th iteration of the Dinkelbach algorithm. Reformulated problem (P4.1) is still non-convex due to nonconvex objective function and constraint (C9) which can be dealt with as follows.
By introducing the slack variables ψ and µ µ µ = {µ[n]} N n=1 such that we can relax the problem (P4.1) to the one with the approximately equivalent but enjoying concave objective function as , ∀n (36e) wherein ω , ∀n} serving as a global convex upper-bound of (11), defined as Note that constraint (36e) must be met with equality at the optimal point, because µ[n] can be otherwise decreased, resulting in an increase of the value of the objective function, which of course, violates the optimality. Plus, we also point out that the objective function, the constraints C10 − C14, and (36c) are now convex. However, the problem (P4.2) is still unsolvable due to the generated extra non-convex constraints (36e) and (36f). Note that the LHS expression of (36e); i.e., summation of norm-square components, is jointly convex w.r.t the variables µ[n] and v[n]. Owing to the fact that the right-hand-side (RHS) of (36e) is convex, since the second derivative of the inverse-square function 1 µ 2 [n] is non-negative; therefore, by replacing the LHS with the corresponding global concave lowerbound using first-order Taylor expansion at the local given point (µ n ) with superscript m indicating the iteration index of fractional Dinkelbach programming, we can reach the approximate convex constraint, associated with (36e), as , ∀n  Note that all the inequality constraints (39d), (39e), and (39f) must also be met with equality at the optimal point, otherwise the optimality is violated. Following the high-SNR approximation, we set ≈ 0 in the subsequent sections for the ease of expositions. We remark the fruitful lemma below. We have the following tight inequalities Proof. Please see Appendix B.
By introducing the slack variables u = {u k [n], ∀k, n}, and using Lemma 3, we can approximate the non-convex problem (P4.3) with a more tractable reformulation given as k,n , ∀k, n} are the value set of slack variables (r, w, u) in the m-th iteration of Dinkelbach algorithm. Finally, since the last constraint is non-convex, we apply [11,Lemma 3] to approximate it with the corresponding convex constraint using the SCA approach, and obtain an approximate convex reformulation of (P4.4) as Result: q , v Initialize feasible point (q (0) , v (0) ) and slack variables, set iteration index m = 0, then It is worth noting that to solve subproblem (P4.5), we have (3N (K + 2) + 1) optimization variables and (3N K + 7N + K + 1) convex constraints. Assuming the accuracy of SCA algorithm for solving this problem is ε 4 , the complexity of solving approximated subproblem (P4.5) for given λ (m) can, therefore, be obtained as O (3N (K + 2) + 1) 2 (3N K + 7N + K + 1) 1.5 log 2 ( 1 ε4 ) . Remark 3. Note that constraints given by (39e) and (39f), being in the form of a x−x 0 2 exp(b x−x 0 ) ≥ y, plus, the expression E = ln(1 + cx −1 + dy −1 ) used in (44c) are proved to be convex; however, they indeed violate the DCP rule-set of the CVX, and so cannot be applied in the optimization model. The former can be handled by rewriting it as 1 , (46) And the latter can be dealt with properly by replacing E-form function appeared in (44c) with t 5 and adding the constraints wherein t 1 −t 5 are some non-zero slack variables, and the logsum-exp function, which is a CVX-approved convex function, defined as LSE(

III-E Overall algorithms and complexity discussion
Having obtained an efficient optimization model for each sub-problem in the previous section, we are now ready to propose iterative algorithms based on sequential block optimization and maximum improvement (MI) or the socalled greedy optimization introduced in [40], summarized in Algorithm 2 and Algorithm 3, respectively. The former is simpler to implement and requires less computations at each iteration. The latter converges faster thanks to a large step-size at each iteration and implementation via parallel computation capability; otherwise, it maybe too expensive. It can be mathematically proved that both algorithms are guaranteed to converge to at least a suboptimal solution. ; l ← l + 1; , run Algorithm 1 with q (l) and v (l) , updating q (l+1) ← q and v (l+1) ← v ; l ← l + 1; 7: Until fractional increase of objective function in (18) gets below the threshold 1 ; 8: Return: (Q opt , P opt , ζ ζ ζ opt )← Q (l) , P (l) , ζ ζ ζ (l) ; , or (q (l+1) , v (l+1) ) whose maximum improvement of objective function given in (18) gets the highest, and keep the remained blocks unchanged; 5: l ← l + 1; 6: Until fractional increase of objective function in (18) gets below the threshold 1 ; 7: Return: (Q opt , P opt , ζ ζ ζ opt )← Q (l) , P (l) , ζ ζ ζ (l) ; Since the feasible solution set of (P) is compact and its objective value is non-decreasing over iteration index l (a similar explanation also applies to the inner Dinkelbach algorithm over the iteration index m), and that the optimal value of minimum SEE is upper bounded by a finite value from the communications engineering perspective [11]. In terms of computational complexity, given L and M be the maximum convergence iteration of the outer overall BCD-SCA algorithm and the inner fractional sub-algorithm, Algorithms 2 and 3 have the complexity of approximately O L(N K) 3.5 log 2 ( 1 ε1 ) + M log 2 ( 1 ε4 ) +LN 3.5 log 2 ( 1 ε2ε3 ) and O L(N K) 3.5 max log 2 ( 1 ε1 ), M log 2 ( 1 ε4 ) . Both are in polynomial time order and applicable to the UAV scenarios.
IV NUMERICAL RESULTS AND DISCUSSION In this section, we provide some numerical simulations to evaluate the secrecy performance of the proposed THz-UUR scheme, and demonstrate the effectiveness of our proposed design in comparison with some benchmarks. Unless otherwise stated, all simulation parameters are given in Table I. Since the initial feasible point is important to use the proposed BCD-SCA-Dinkelbach based algorithms and significantly impacts their convergence performance, we explain how we can obtain initial feasible UAV's trajectory, velocity, transmission powers, and user scheduling. The initial UAV's trajectory is assumed to be circular centered at the BS's location with radius R u = q b − q I , provided that UAV's instantaneous velocity constraint C12 is satisfied, and T ≥ T min Vmax , where T min cir is the minimum required time for circular trajectory. However, if T min Vmax (i.e., at least cyclic trajectory was possible with minimum required time T min cyc ), then one could use any cyclic shape as long as C10 − C14 are satisfied. Here, we consider a Piriform trajectory with discretized equations given by q i = [x i ; y i ] with y i = A y (1 − sin(t)) cos(t) and x u = R u (sin(t) + 1)/2 in which t 1×N indicates the linearly spaced vector in [ π 2 , 5π 2 ]. Further, the constant A y can be obtained efficiently via a simple 1D search in the range of [R u , 0] or simply set to zero. The UAV's initial velocity vector v i is then followed by Having obtained an initial feasible UAV's trajectory and velocity and the UEs are scheduled equally (e.g., N K times each), i.e., ζ ζ ζ i is obtained such that the constraint C1 holds.
After identifying the initial feasible point for the iterative optimization algorithms, we consider different benchmark schemes, all of which are detailed below and labelled in the following figures, to demonstrate the superiority of our proposed minimum SEE-based optimization algorithms.
• SEE-Seq: minimum Secrecy Energy Efficiency optimization scheme using the Sequential BCD-based subproblem maximization as given in Algorithm 2.   • ASR-Seq: Optimizing minimum ASR given in (15) while ignoring the UUR's flight power limit using the Sequential BCD approach to iteratively improve Q, P, ζ ζ ζ. Fig. 2 depicts the convergence of the proposed iterative algorithms. We can see that both benchmark schemes SEE-FTrj and SEE-FPow converge quickly; however, they can only achieve significantly lower SEE performance than the proposed trajectory using joint design of power control and user scheduling. Specifically, SEE-MI converges relatively faster than SEE-Seq, i.e., 13 against 28 iterations, at the cost of slightly lower minimum SEE than that of its counterpart. However, they both achieve approximately 68% SEE improve-    Fig. 3 illustrates UUR's trajectories using different optimization algorithms. We see that the optimized trajectories are much more complicated than the initial circular one with the counterclockwise direction. Notice that UUR should fly towards UEs' locations to obtain data with low power. This, in turn, can significantly increase the chance of information leakage due to a stronger wiretap link and less effective BS's jamming. Thus, UUR prefers to stay not too far from the BS. Overall, we see that the path planning makes UUR adjust trajectory through the best possible path, efficiently forming the distances between the UUR, selected UEs, and the BS such that the trade-off in the channel conditions for the friendly jamming transmission in the first phase as well as the aerial relaying in the second phase of transmission improve the secrecy performance. Further, we observe that the SEE-based trajectories are smoother than that of the ASR-Seq scheme, implying possibly a lower flight power consumption of UUR. The SEE optimization demands this in contrast to the ASR-Seq design where the UUR's velocity might harshly fluctuate for the minimum ASR (mASR) improvement if required. We also note that when the initial circular trajectory is impossible due to significantly low mission time, e.g., T = 5s, and owing   to the UAV's physical system limitations, the crucial task of path-planning can be efficiently designed as shown in Fig. 4. It should be mentioned that the curve belonging to the "SEE-FTrj" does represent the initial feasible cyclic trajectory based on the Piriform, and the other curves illustrate the optimized UUR's trajectory according to the different algorithms. Fig. 5 illustrates mASR and the average flight power consumption (AFPC) against iteration indices for different schemes. It is crystal clear that for the SEE-based algorithms, the mASR and the AFPC performances tend to be nondecreasing and non-increasing, respectively. In contrast, for ASR-Seq scheme, the AFPC first decreases and then increases until convergence after 30 iterations. We also note that this scheme can achieve slightly higher mASR performance than our proposed schemes but at the cost of significantly lower SEE (43.13 Mbits/Joule). Fig. 6 is plotted to demonstrate how the UAV's velocity (Vel.) and the instantaneous flight power consumption (IFPC) can be adjusted over time for SEE improvement using different algorithms. We observe that all SEE algorithms, except "SEE-FTrj", make UUR fly with roughly less speed variation for a relatively more extended period of time (e.g., from 3s to 8s) to satisfy mission requirements as well as improve the SEE performance. However, due to having complicated function of IFPC w.r.t the UAV's velocity given in (11), UUR starts at a high initial speed to fast reach the targeted location, but not at maximum speed or hovering for the sake of efficient power consumption purposes. Fig. 7 illustrates the joint power allocation and user scheduling vs. time for different algorithms. The sub-figure 7d represents the non-optimal but feasible power allocations and user scheduling adopted for initialization of all the algorithms. Initially, UUR is very close to UE 5 but far from the BS. Hence, UE 5 is scheduled due to a possibly better channel condition than the others, and the BS jams in high power while UE 5 keeps low power. For SEE-FTrj, UUR follows the circular trajectory while maintaining the same distance from the BS that has a constant jamming power. In contrast, subfigures 7a, 7b, and 7e show that at initial stage, UE 5 increases power when UUR heads towards the BS and the BS reduces jamming power. Further, these UEs are scheduled unequally, but during their scheduling, except UE 5 , they need to utilize their maximum transmission powers for sending information, and the relaying power slightly fluctuates around p ave u . Finally, Fig. 8 depicts how the SEE performance varies when the molecular absorption coefficient of THz links changes from a f = 0.005 to a f = 0.025, which is proportional to the carrier frequency. When the mission time increases from T = 5s to T = 13s, the SEE performance improves due to more time for secure communications and adjusting flight parameters. It also demonstrates that the larger the molecular absorption coefficient, the lower the SEE performance for low mission time (T = 5, 8, 9s) due to higher propagation loss arising from severe molecular absorption. However, it is worth pointing out that the increased propagation loss results in the reduction of not only UUR's information leakage, but also BS's reception quality. The overall trade-ff between these two phenomena, therefore, results in the fact that the SEE performance does not get monotonically decreased as a f increases, according to curves T = {11s, 13s}.

V CONCLUSIONS
In this paper, we investigated the challenging task of designing an energy-efficient THz-UUR system for secure and periodically data delivering from multiple ground UEs towards the BS. For the fairness of QoS amongst the UEs, a minimum SEE maximization problem was formulated, by which the fundamental system parameters are designed to improve the overall system secrecy and energy-efficiency performance. This was formally posed as a challenging mixed-integer nonconvex nonlinear maximin optimization problem. We then embarked on tackling the nonconvexity of the formulated problem and proposed low-complex BCD-SCA-Dinkelbach based iterative algorithms to solve it suboptimally with guaranteed convergence. Simulation results confirmed the fast convergence of our proposed algorithms, demonstrated significant SEE performance improvement than the other benchmarks, and provided insightful results in the optimized system parameters such as UUR's trajectory and velocity pattern as well as communication resource allocations, including transmit power profiles and user scheduling. Also, the effects of mission time, and molecular absorption factors arising from the THz links on    the system SEE performance have been examined. As future work, we will deeply investigate the dynamic topology of aerial platforms with more practical THz channel modeling for intelligent UUR systems.
APPENDIX B PROOF OF LEMMA 3 Computing the gradients of given functions w.r.t x and y yields We can verify that the first-order and second-order determinants of H 41 and H 42 are all non-negative, and therefore, the Hessian matrices are positive semi-definite (H 41(2) 0), indicating that functions f 41 (x, y), f 42 (x, y) are jointly convex w.r.t x and y. Further, the convexity of f 43 (x, r) and f 44 (x, p) follows from (B.6) and [15,Lemma 1]. Given these functions are all convex, one can use the first-order Taylor expansions at points x 0 and y 0 to reach the global tight lower-bounds and inequalities in Lemma 2. The proof is completed.