{'fraction_non_alphanumeric': 0.0365567307899442, 'fraction_numerical': 0.009133797203851705, 'mean_word_length': 4.220060721916114, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Green Tethered UAVs for EMF-Aware Cellular Networks

A prevalent theory circulating among the non-scientific community is that the intensive deployment of base stations over the territory significantly increases the level of electromagnetic field (EMF) exposure and affects population health. To alleviate this concern, in this work, we propose a network architecture that introduces tethered unmanned aerial vehicles (TUAVs) carrying green antennas to minimize the EMF exposure while guaranteeing a high data rate for users. In particular, each TUAV can attach itself to one of the possible ground stations at the top of some buildings. The location of the TUAVs, transmit power of user equipment and association policy are optimized to minimize the EMF exposure. Unfortunately, the problem turns out to be mixed-integer non-linear programming (MINLP), which is non-deterministic polynomial-time (NP) hard. We propose an efficient low-complexity algorithm composed of three submodules. Firstly, we propose an algorithm based on the greedy principle to determine the optimal association matrix between the users and base stations. Then, we offer two approaches, a modified K-mean and shrink and realign (SR) process, to associate each TUAV with a ground station. Finally, we put forward two algorithms based on the golden search and SR process to adjust the TUAV's position within the hovering area over the building. After that, we consider the dual problem that maximizes the sum rate while keeping the exposure below a predefined value, such as the level enforced by the regulation. Next, we perform extensive simulations to show the effectiveness of the proposed TUAVs to reduce the exposure compared to various architectures. Eventually, we show that TUAVs with green antennas can effectively mitigate the EMF exposure by more than 20% compared to fixed green small cells while achieving a higher data rate.

it suffers from high path loss. Recently, reconfigurable intelligent surfaces (RISs), which can reflect the incident electromagnetic wave toward the specified direction, have been proposed to reduce exposure [12]. However, in order to operate on a large scale, several RISs need to be installed and optimized jointly, which can be a challenging task. Moreover, the signal processing capabilities of the RIS are extremely limited.
Another candidate for widespread deployment of 5G and beyond communication systems is small cell offloading, which has proven to improve the performance and enhance the energy efficiency [13], [14]. The main idea of this method is to unload the user traffic from the macro cell to small cells (SCs). Thoroughly, a significant reduction in the UL radiated power can be achieved [14]. The radiation can be further reduced by considering the decoupling between UL and DL [15], where green antennas, i.e., receiving only radio, are used to assist the UL [16].
Although the use of fixed deployed SCs and RISs can effectively reduce EMF exposure, the network cannot adapt to user distribution changes [12], [16]. Therefore, a more flexible deployment with airborne-SCs can assist the existing cellular infrastructure to provide users with better service quality, e.g., coverage and capacity [17]- [21].
In [22], the authors proposed an efficient cell-based allocation method, which provides an optimized unmanned aerial vehicle (UAV) positioning to enhance the performance of communication systems. Moreover, UAVs can hover over geographical regions that experience heavy traffic conditions due to natural disasters or mass events to assist the existing infrastructure in providing users with better service quality [23]. But the limited onboard energy and flight time impose a critical challenge for the deployment of UAVs [24]. In this regard, TUAVs can be a viable alternative, as they can be powered and backhauled through cables connecting them to ground stations (GSs) usually located on the rooftop. TUAVs can outperform regular UAVs in terms of coverage and capacity [25], [26]. Nevertheless, the exposure in the DL direction can increase. Also, the position of the TUAVs and the resource allocation problem on a large scale have not been investigated.
In contrast to existing works that mainly consider analyzing the coverage, our goal is to design EMF-aware cellular networks with low complexity algorithms to minimize the exposure and guarantee a target QoS. To the best of our knowledge, EMF-efficient architectures and association algorithms involving TUAV carrying green antennas have not been discussed in the literature.

B. Contributions
In order to alleviate some of the issues in the common adopted architectures with UAVs, e.g., the reduced mobility, limited power capacity, excessive radiation in DL, and backhauling, we propose a system with green TUAVs to minimize the EMF exposure by optimizing the transmit power, users' allocation, and the location of the TUAVs.
More precisely, we propose an architecture where there are several GSs that TUAVs can attach to. Moreover, each TUAV has a limited number of RBs. This problem tends to be a mixed integer non-linear programming (MINLP); hence, we design several efficient algorithms to solve the optimization problem. The main contributions of the paper can be summarized as follows: • We propose a novel architecture for cellular networks that exploit green TUAVs as an efficient tool to minimize the EMF exposure. • We design alternate optimization algorithms to optimize the users' association, deployment of TUAVs, and fine positioning of the TUAVs within their hovering area. The complexities of the proposed algorithms are analyzed.
• We consider the dual problem, i.e., maximizing the average UL rate with a constraint on the exposure, e.g., the ICNIRP limit. Then, we present algorithms for the efficient deployment of the TUAVs for the dual problem.
• We assess the performance of the proposed scheme against various architectures, e.g., different schemes with TUAVs and fixed SCs that assist the UL and/or DL. In the lower right corner, the model of green TUAV and BS are shown separately.

C. Notations and Organization
The main notations used in the paper are listed in Table I. The remainder of the paper is organized as follows : Section III explains the considered system model. The problem formulation for minimizing the exposure with achieving a target QoS is given in Section IV. Section V details three submodularity of the objective functions and heuristic algorithms. Next, we present the dual problem, where we propose a novel design of cellular networks with an EMF constraint in Section VI. Then, simulation results are provided in Section VII. Lastly, Section VIII concludes the paper with a few remarks.

III. SYSTEM MODEL
We consider a TUAVs-enabled system consisting of one macro base station (BS), K users and M TUAVs with green antennas, i.e., receiving only antennas [16], as show in Fig. 1. We also consider N GS placed in fixed locations at the top of some buddings and connected to links to provide backhaul and power to the TUAVs, where each TUAV can choose one GS to connect to. We consider that the BS is located in the center of the area and is responsible for resource management. 1

A. Channel Model
In the considered scenario, each user can receive two main signals, i.e., LoS signal and strongly reflected NLoS signal from gNBs [29], the corresponding path loss for the LoS and NLoS propagation can be written as respectively, where f c is the carrier frequency, c is the speed of light, r kj is the distance between the gNB j and the user k, α l and α n are the path loss exponent for the LoS and NLoS links, respectively. The symbols F l and F n denotes the random variables (RVs) representing the power loss due to large scale fading and shadowing with mean η l and η n [30]- [32]. Moreover, G l and G n are RVs accounting for the power loss due to small scale fading with mean unity [33].
We adopt the channel model based on the probabilistic LoS model, averaged over the large and small scale fading [30]- [32]. The LoS and NLoS links have different probabilities of occurrence, where the probability of LoS can be written as where d kj is the distance from the projected position of gNB j to user k, and a and b are constants that depend on the environment. For the probability of the NLoS signal, we have p NLoS kj = 1 − p LoS kj . Therefore, the average path loss can be writte as

B. UL Association
In the system, TUAVs with green antennas are designed to reduce exposure to EMF. In the DL, all users are connected to the BS. For UL, some users are connected to TUAVs, while the rest can be associated with the BS.
Furthermore, we introduce a binary variable ǫ kj , which is equal to 1 if user k is associated with gNB j and 0 otherwise, i.e., We also consider that the number of available spectrum resources at gNB j is W j,max RBs, and each RB has a bandwidth of B that can be allocated to a single user for either UL or DL communications, i.e., the maximum number of users that can be associated with that TUAV is W j,max . Hence, the following conditions should be satisfied K k=1 ǫ kj ≤ W j,max , ∀j ∈ J {0, 1, · · · , M }, For convenience, we consider that the index j = 0 refers to the BS, while the other elements represent the TUAVs.

C. Deployment of TUAV
For the deployment of TUAV, each TUAV can only be connected with one GS at a certain time. More precisely, where ϑ mn is a binary variable that equals 1 only if the m th TUAV is attached to the n th GS.

IV. PROBLEM FORMULATION
In this section, we describe a metric usually adopted to quantify the EMF exposure and propose a power allocation strategy for the devices. Then, we formalize the optimization problem such that the EMF exposure due to UL transmissions is minimized while guaranteeing a certain QoS in terms of the users' UL rate requirement.

A. EMF Assessment
To evaluate the EMF radiation in a cellular network, we adopt an exposure metric called EI [34], which is able to model the population exposure to EMF, accounting for various elements such as different technologies, environments, and usages. The EI can be written as where the exposure in each case can be represented as a function f of the mean emitted powerP TX from the device, and mean received power densityS RX , and the SAR reference levels SAR UL and SAR DL . The EI depends on the population age group l (e.g., adults, seniors, and children), environments e (e.g., indoor and outdoor), radio access technologys (RATs) and layers r (e.g., WiFi, 5G, narrowband -Internet of things (NB-IoT), and macro, micro), usage types g (e.g., voice, data, and machine-type communication).
In this paper, we consider one RAT which is next-generation cellular network and two layers: i) macro BS; ii) TUAV carrying a SC with green antennas. Then for simplicity, we consider adults as the population, while data and voice are the considered usage types. More precisely, now the EI for both UL and DL can be given by with EI UL representing the UL exposure given as where SAR UL k is the whole body SAR or localized SAR normalized to unit transmit power, which depend on the required service and the posture (how the device is hold), P UL k is the transmit power of users k. For the DL exposure EI DL , it can be obtained as where SAR DL is the DL SAR value which is normalized to a reference received power density of 1 W/m 2 and K is the number residents in the region (including users and non-users), i.e., K ⊆K, P DL jk is the received power density at resident k from gNB j.

B. Transmit Rate and Power
Let P kj be the transmit power of user k when connected to the j th gNB, his UL rate can be expressed as where σ 2 = k B T B is the noise variance with k B being the Boltzmann constant and T being the temperature in Kelvin. 2 Next, we focus on DL and assume that gNBs are powerful enough to meet DL rate requirement of the connected user, and total transmit power of gNB j can be given as where the subset K j represents users connected to gNB j. The required DL power of each user is a function of the DL rate R DL kj and the path loss L kj .

C. Complete Optimization Model
In (9), we consider a general model where both UL and DL exposure are accounted for. The total EMF exposure is mainly dominated by UL due to the proximity of the device to the body [3], [14], as shown also in Section VII.
Therefore, we only focus on the UL exposure, and our objective is to optimally deploy the TUAVs and associate the users in order to minimize the UL EI. Moreover, there should be a minimal power received by gNB to satisfy the required rate, which depends on the path loss and transmit power. In fact, the transmit power should compensate for the path loss, which depends on the association between TUAVs and GSs ϑ mn , the association between users and gNBs ǫ kj , as indicated by (12). Also, the path loss depends the distance between the TUAV and the users, r kj , as in (4), which, by its turn on, depends on the tether length T m , elevation angle θ m , and azimuth angle ϕ m of the tether, as shown in Fig. 2.
Let γ {ϑ mn , ǫ kj , T m , θ m , ϕ m } , j ∈ J , k ∈ K, m ∈ M, n ∈ N represents all optimization variables, i.e., the association between users and gNBs, association between TUAVs and GSs, tether length, and inclination angles. Then, our optimization problem can be formulated as where the objective function is the exposure index that can be seen as a weighted sum of the users' transmit power such that the weight of each user is his constant reference SAR, SAR UL k . For the constraints, (14b) guarantees achieving a minimum required UL data rate, (14c) limits the transmit power of users due to the hardware restrictions, and (14g)-(14i) ensure that TUAVs are within their hovering area with θ m > θ min for safety reasons. The problem P0 is mathematically challenging because it includes integer and continuous variables, as well as non-convex objective function and non-convex constraints, which makes it a MINLP.

V. OPTIMIZATION ALGORITHMS
In order to solve the target optimization P0, we propose an alternate optimization algorithm, where we split the problem P0 into three sub-optimization programs, (i) finding the association matrix, (ii) determining the deployment of TUAVs, and (iii) adjusting the tether length and angles for each TUAV. However, it is still complicated to solve each step due to the numerous combinations or non-convex continuous problem. Thus, we propose low complexity heuristic algorithms to find the optimal 3D placement of TUAVs and users association to minimize the EMF exposure.

A. Power Allocation
The actual rate should be larger than or equal to the required rate according to constraint (14b). Since the EI monotonically increases with the users' transmit power, the optimal power is the one that achieves the required rate with equality. However, there is also a limitation in power due to hardware limitations, as seen in constraint (14c). Therefore, the transmit power P k (γ) can be allocated as where P max is the maximum transmit power of the users, and j * k is the index of serving gNB for user k, i.e., ǫ kj * k = 1. It can be seen that power allocation scheme consider the minimum between required and maximum transmit power, where the power constraint is more dominant than that of rate, because it is a hardware constraint that the mobile equipment cannot exceed. If the required transmit power to achieve the rate is larger than the maximum transmit power, the problem will be infeasible. Therefore, the achieved rate will be less than the required one.

B. Optimize the Association between Users and gNBs
For arbitrary given locations of the TUAVs, we focus on how to associate each user with the appropriate gNB for the UL. Under this conditions, the association problem can be formulated as where we consider the power allocation strategy in (15). The problem is also non-deterministic polynomial-time (NP)-hard, as we have to go through all the combinations between users and gNBs. In order to solve the problem, we propose an algorithm, which starts with each user connected to the best gNB, i.e., the one corresponding to the minimum path loss. Then, for the gNBs with an excessive number of connected users, we will change the connection of some users. That is, some users will be served by other gNBs which are not fully occupied, rather than the best one. Moreover, the users to be moved are selected such that the induced increase in radiation, which results from the non-optimal selection of the gNB, is minimal.
More precisely, the increase in exposure when user k is connected to the j th gNB instead of the best choice j * can be quantified asḡ kj e kj − e kj * , where e kj SAR UL k P kj is the EMF exposure generated by user k when assigned to gNB j. Let us defineK as the set of users connected to excessive occupied gNBs, andJ as the set of gNBs that are still not fully occupied. If the number of connected users exceeds the resource limit at some gNBs, we search for the user k c ∈K with the smallest induced increase in radiation (i.e., smallestḡ kj for all k ∈K and j ∈J ) to be moved first. This user k c will be moved to the gNB that minimizes the induced radiation when the user moves. This process is repeated until each gNB is associated with a number of users that do not exceed the available resources.
For the next stage, where we optimize the location of the TUAVs, we need to quantify the role that each TUAV contributes for reducing the EMF. This contribution of the j th TUAV to reduce the radiation can be defined as g j k∈Kj g kj , ∀j ∈ M, where g kj e k0 − e kj is the gain, when user k is connected to the j th TUAV instead of the BS, as described in Algorithm 1.
Next, we analyze the complexity of Algorithm 1, which is mainly composed of two parts: i) the initial association, ii) users' re-association due to the limited resources at TUAVs. The algorithm's complexity mainly comes from the first stage, that is in the order of J × K. In the second stage, from the users' perspective, each user only needs to be changing his association once at most, and the number of available gNBs is also less than J. Hence, the complexity of Algorithm 1 is O(J × K).

C. Optimize the Deployment of TUAVs to GSs
Unlike regular UAVs, the freedom of TUAVs is limited by the predefined locations of the ground stations, tether length and inclination angles. In order to reduce the exposure, a proper location for the TUAVs should be determined, i.e., TUAVs' association with GSs. The deployment problem is a MINLP, as we need to consider all possible configurations between TUAVs and the GS as well as the inclination angle and tether length of each TUAV. For given tethered length and inclination angles, e.g., TUAVs are located in the center of their hovering areas above the GS, the TUAV association to GSs problem can be written as The Problem P2 is a binary integer program for variables ϑ mn and ǫ kj . One approach to avoid exhaustive search over all possible configurations is to select the locations that are most likely to provide satisfactory performance, Algorithm 1 Association matrix 1: Input W j,max , j ∈ J , location of gNB j ∈ J , location of users k ∈ K 2: W j = 0, ∀j ∈ J , generate initial number of associated users to gNB j 3: initial association matrix ǫ kj = 0 for each k ∈ K, j ∈ J 4: compute the quantified EMF exposure and gain, e kj and g kj for each k ∈ K, j ∈ J , respectively 5: for k = 1 : K do {Initial assignment of users } 6: find j * k s.t. j * k = arg min j∈J (e kj ) (i.e., find the best serving gNB for user k, regardless of resource constraint) 7: , move one user to a non-fully occupied gNB, the corresponding increment of this alteration is minimal.) 12: W jc = W jc + 1, and ǫ kcjc = 1 13: W j = W j − 1, and ǫ kcj = 0 14: end while 15: end for 16: calculate TUAV's contribute to reducing radiation g j = k∈Kj g kj , ∀j ∈ M 17: Output e kj , ǫ kj , g j , ∀k ∈ K, ∀j ∈ J such as the barycenter of users' concentrations weighted by their exposure. 3 Another method is a divide-and-conquer like algorithm, where we iteratively prune the search space on the premise of ensuring that each iteration is moving towards a better result. For this purpose, we propose two heuristic algorithms, based on SAR-aware K-mean and shrink and realign (SR) process [37], respectively. Both of them are composed of three steps: (i) deploying randomly the TUAVs, (ii) changing the position of TUAVs, (iii) searching for the nearest GS. In the following, we describe the two proposed algorithms in detail.
1) Modified K-mean Algorithm: Initially, TUAVs are randomly scattered over the entire area, while the location of the BS remained unchanged in the middle of the area. Then, the users are associated to the gNBs according to Algorithm 1 and the gain of TUAVs at the initial iteration, g i j , j ∈ M, of their deployment is also computed. Let us define x mid j , y mid j as the SAR-weighted barycenter of users associated with the j th TUAV. We move the UAVs to the barycenters and run Algorithm 1 again to obtain the gains of TUAVs if they are moved to the new positions, g mid j , j ∈ M. Then, we pick the coordinates based on the value of gains, i.e., select either the current location of the TUAV or the user's weighted barycenter depending on the location that can reduce more the exposure. This process is repeated until the difference between the i th and the (i + 1) th coordinates are less than the tolerance δ, or reaching the maximum number of iteration I max . Finally, each TUAV is connected to nearest GS from unoccupied GSsN , while the location of the BS remains the same, as described in Algorithm 2.
Algorithm 2 Determine TUAVs' deployment by K-mean 1: Input W j,max , j ∈ J , location of users k ∈ K, I max 2: i = 1 3: initial connection matrix ϑ mn = 0 for each m ∈ M, n ∈ N 4: generate coordinates x i j , y i j that TUAVs are randomly scattered in the 2D space, and the BS is fixed in the center. Run Algorithm 1 with gNBs' coordinates x i j , y i j , j ∈ J to get EMF exposure e kj and association matrix ǫ kj for users and gain g i j for TUAVs 7: x mid j = k∈K j x kj e kj k∈K j e kj , j ∈ M Run Algorithm 1 with the coordinates of gNBs x mid j , y mid j , j ∈ J to get gain g mid j for TUAVs 10: 2) Algorithm Based on 2D SR Process: Initially, TUAVs are uniformly scattered over the entire area, while the location of the BS remained unchanged in the middle of the area, q i j , j ∈ J . Then, we generate initial next position candidates q i j,t , ∀t ∈ T {1, 2, · · · T } as a circle with radius r i around all current locations of the TUAVs to form initial positions. Next, we will select one TUAV and iterate over the candidate locations of it, solve the associations by Algorithm 1, and calculate the total EMF exposure EI SR (q i j,t ) k∈K j∈J e kj ǫ kj , while keeping all the other TUAVs fixed. Then, we find the best next TUAV candidate and compare it to where the TUAV is located to decide whether to move it or not, i.e., choose the location where there is less radiation over the area when this TUAV is in that location. Furthermore, we will move TUAVs around one by one in this way, then generate new candidates on a circle of radius r i+1 = r i /2 around each local solution. We repeat this process until the size of sample space decrease below a certain threshold r min . Eventually, we move the TUAVs directly over the corresponding nearest GS gs n , as described in Algorithm 3.

Algorithm 3
Optimize the deployment of TUAV through 2D SR process 1: Input W j,max , j ∈ J , location of users k ∈ K, r min 2: i = 1 3: initial connection matrix ϑ mn = 0 for each m ∈ M, n ∈ N 4: generate initial location of gNBs q i j , j ∈ J that UAVs uniformly scattered and the BS fixed in the center 5: generate initial candidates in a circle of radius r i around each TUAV q i j,t , t ∈ T 6: while r i > r min do 7: for j = 1 : M do 8: find q i j,t0 s.t. q i j,t0 = arg min t∈T (EI SR (q i j,t )) (i.e., run Algorithm 1 to find association matrix ǫ kj , j ∈ J , k ∈ K and calculate the total exposure EI SR ) 9: the process of finding the SAR-weighted barycenter and moving TUAVs to new location is the most complex step. In each iteration, its complexity mainly comes from the operation of Algorithm 1. After that, the number of iterations is finite, its value is much smaller than the number of users. Hence, the complexity is O(J × K). Moreover, for Algorithm 3, the process of gradually narrowing the search radius is the most complex. More specifically, before each reduction in the search radius, we need to search over all candidate locations of each TUAV, during which we run Algorithm 1 constantly to calculate the total radiation. So the complexity is O(J 2 × K).

D. Optimize TUAV's Position within the Hovering Areas
Through the above two-step iteration, we can fix the deployment of TUAVs and association matrix. Now, we optimize the location of TUAVs on a small scale in their hovering zone. The goal is to find optimal inclination angels and tether length of each TUAV, and this sub-optimization can be written as subject to: (14b), (14c), (14g), (14h), (14i).
In order to avoid entanglement with surrounding buildings and ensure safety, the inclination angel θ of the tether has a minimum value θ min , which is related to environments [25]. Then, the distance between TUAV and GS, T m , should be smaller than the maximum tether length T max , and there are no additional restrictions on azimuth ϕ k .
It can be seen that problem P3 includes continuous variables, i.e., θ m , ϕ m and T m , and it involves a non-convex objective function. We propose to solve this sub-optimization via two iterative algorithms, i.e., golden section search and SR process. For both algorithms, we consider alternate optimization, where each TUAV is individually optimized while keeping the users' association unchanged, as computed in preceding algorithms.

1) Algorithm Based on Golden
Search: At first, we consider a random tether length T m for the m th TUAV.
Then, we find the SAR-weighted barycenter of its associated users, i.e., the projection in the horizontal x-y plan.
If that location is not in the hovering area, we will choose the nearest point in the flyable zone. Next, we compute the EI at this height when the TUAV is moved to the selected point O m (x m (h m ), y m (h m ), h m ), i.e., where h G is the height of GS at the top of the building, and T max is the maximum tether length. After that, we could narrow down the radius based on the golden search, where v = is the golden ratio. The initial heights are the upper and lower bounds, i.e., h u = h G + T max and h l = h G , respectively. Also, we will pick two golden Then, we calculate exposure at this two altitude EI(O u ′ ) and EI(O l ′ ). Based on the radiation of this two points, we will narrow down the search area. Then, this process will iterate over and over to narrow the limits until the interval length is smaller than the tolerance height h min , as described in Algorithm 4. We opt for the golden search method, as it converges fast and is less computationally demanding than an exhaustive search over a large number of position candidates [38].

2) Algorithm Based on 3D SR Process:
The proposed algorithm is based on SR process to optimize each TUAV location individually by choosing the proper location among T ′ candidates around the current location of the TUAV, denoted by q i m , at iteration i. Let us define q i m,t as the t th potential position for the location of the m th TUAV on polyhedron around q i m with radius r i . Then, we calculate radiation EI q (q i m,t ) k∈Km e km when TUAV at those candidate locations. Then, we select the location, among the candidates, that minimizes the radiation EI q . After that, by SR process, we generating T ′ new candidates on a polyhedron with radius r i+1 = r i /2 around the TUAV.
We repeat this process until the radius is less than the precision r ′ min , as described in Algorithm 5.

E. Overall Complexity
Generally speaking, there are two options for determining the deployment of TUAV, modified K-mean and algorithm based on 2D SR process, and the complexity of them is O(J ×K) and O(J 2 ×K), respectively. Therefore, only considering the complexity of algorithms, K-mean is a better choice at this step. After the deployment of   general, the number of users K is usually much larger than other parameters, so the maximum complexity of the entire process should be determined mainly by the number of users. In this regard, we propose the use of K-mean for connection between TUAVs and GSs, and SR process to adjust TUAVs position in a small range, as they have the lowest complexity and best performance in terms of minimizing radiation as shown in numerical results.

VI. DUAL PROBLEM: CELLULAR SYSTEM DESIGN WITH EMF CONSTRAINTS
In previous sections, we consider the problem of minimizing the exposure, while achieving a target rate. Another architecture of interest for cellular operators is the dual problem, where we aim to maximize of the rate under a constraint on the users' exposure to EMF. More precisely, our objective is to optimally deploy TUAVs in the space find q i m,t0 s.t. q i m,t0 = arg min t∈T ′ EI q (q i m,t0 ) (i.e., choose the best candidate position which is in hovering area) 8: where SAR limit is the whole body or local SAR threshold. The constraint (20b) guarantees that the EMF exposure of each user is less than a pre-specified limit, e.g., the ICNIRP limit, to avoid the health risks, while other restrictions remain the same. In the following, we propose an algorithm to solve the problem.
In order to improve the sum-rate, the user equipment (UE) can send the maximum allowable power that satisfy constraint (14c). Nevertheless, the induced exposure to EMF can be escalated above the required level, violating the constraint (20b). Hence, an additional limit on transmit power is added to satisfy both of the constraints, i.e., where transmit power P k (γ) is the minimum between the maximum transmit power P max and the power to ensure the SAR limitation SAR limit /SAR UL k . In fact, each user k can have his own reference SAR, SAR UL k , depending on several factors such as postures (standing or sitting), usages (voice or data), and device type (wearable, handheld, or in-body).
Although the objective function in P4 is different from the original one P0, the dual problem is also a MINLP.
Therefore, the proposed algorithms in section V are still applicable and we only need to make some minor 18 adjustments as follows: • We consider a new power allocation strategy as in (21) rather than (15).
• For Algorithm 1, the radiation e kj , is replaced with the minus of the rate −R kj , and the gain is redefined as • For Algorithm 2, we change the function from SAR-weighted to Rate-weighted, i.e., using rate −R kj to replace e kj in the formula for calculating weighted barycenters.

VII. NUMERICAL RESULTS
In this section, we present some ii) Poisson cluster process (PCP), were each cluster represents one hotspot with a radius 100 m, with average number of usersK/6 , and the average number of clusters in the area is four. In simulation, all the metrics used are averaged over 1000 iterations. We consider that active users fall within two categories according to their usage, i.e., voice and data users. Each category has distinct UL rate requirements and SAR reference for the EMF exposure. In fact, since the mobile phone is closer to the person when making a call through the handset directly, the induced whole-body SAR of voice usage is larger than that of data usage. For instance, we consider a reference SAR, SAR UL k ∈ {SAR v , SAR d } with SAR v = 0.0047 and SAR d = 0.0037 W/kg per unit transmit power in the sub-5 GHz band for voice and data users, respectively [34, Table 27].
For convenience, the simulation parameters are summarized in Table II, unless specified otherwise. In the following figures, for the scenario named BS only, we consider a single BS in the center of the area, and all the users are connected to it for both UL and DL.

A. Performance of the Proposed Algorithms on Small Scale
We move our attention to the performance of the proposed algorithms for each sub-optimization process. In section V, it can be seen that we solved the optimization problem in three steps. Firstly, we determine the association between the user and gNBs for the given position of TUAVs, then ascertain the deployment of TUAVs in a 2D plane, and finally, slightly adjust the tether length and angle of each TUAV. In the process of comparison, we consider the results of brute force search as a baseline. Therefore, in this simulation, we only consider two fixed TUAVs evenly distributed, each of which has two RBs. In addition, there are a BS in the center of the region,  N = 9 evenly distributed GSs, andK = 60 residents. The number of active users K ranges from 4 to 20, and all active users are data users whose required UL rate is 50 Mbps.

1) Optimize the Association between Users and gNBs:
In the beginning, we compare the effectiveness of different methods for determining the association matrix. In addition to the results of our proposed algorithm, we also show those of BS only, randomly associating the user with gNBs, and brute force approach. As can be seen from Fig. 3, randomly associating users with gNBs can even be counterproductive. The result of our proposed algorithm almost overlaps with that of brute force. Therefore, the association policy in the following simulations is determined via our proposed algorithm. 2) Optimize the Deployment of TUAVs to GSs : Fig. 4 reports the numerical evaluation of radiation, considering associating the TUAVs with proper GSs. We compare the results of BS only, random deployment of TUAVs, and traversing all combinations (brute force search) with the results of the two proposed algorithms. As shown in Fig. 4, the results of the modified K-mean and SR algorithms almost coincide, and there is not much difference between them and that of brute force. Furthermore, these two algorithms can provide an additional 15% gain compared to the TUAVs' random deployment. Combined with the previous complexity analysis, we chose the modified K-mean with less complexity to determine the TUAVs' deployment in the following results.

3) Optimize TUAV's Position within the Hovering:
In the final step, we consider each TUAV separately and only move them in a small area around the GS. Furthermore, we also compare the results of the two proposed algorithms with the other three baselines, BS, random selection, and brute force search over a gridded hovering zone. As can be seen from Fig. 5, the optimization over the fine location of the TUAVs in the hovering area can result in a gain, albeit not as large as the previous two. Moreover, although the results of the two methods are relatively close, the gain brought by the adopted golden search is slightly smaller than that of the 3D SR process. So the latter approach will be used in subsequent simulations.
Another alternate optimization iterations can be considered after solving the problem by three-step sub-algorithms.
In this iteration, the optimized variables of the last iteration are provided as the initial values of the optimization variables. Such alternate optimization should further enhance the performance, at the expense of additional complexity.

B. UL and DL Exposure with TUAVs
In the following, we analyze the effects of different types of TUAVs (i.e., various UL and DL decoupling methods) on reducing EMF exposure. In each scenario, TUAVs have N = 25 GSs to choose from. We haveK = 120 residents in the area, K out of them are randomly selected as active users, who have two different ways of using the phone, voice (20% of users) or data (80%). We consider that the maximum available RBs for each TUAV is W m,max = 6, while BS has sufficient resources to allocate all users in the area. If TUAVs are used to densify(i.e., assist) both UL and DL, then each user needs to occupy two RBs, and the maximum number of users connected to each TUAV is W m,max /2. There are four scenarios, of which BS only is previously described, and the other three involve M = 4 TUAVs, and they are described as follows: • Green TUAVs (Assisting UL): This scenario is our proposed network architecture, where the location of BS is fixed like last scene, and we additionally add green TUAVs which can only receive information from users and forward it to BS through cable. Moreover, all users receive signals from the BS, and in UL up to 6 users are connected to each TUAV. Therefore, we run Algorithm 1 to adjust association matrix, while operate Algorithm 2 to regulate the position of TUAVs in 2D space. After the above two steps converge, we will adjust the tether length and angles of each TUAV through Algorithm 5.
• Regular TUAVs (Assisting UL and DL): In this scenario, the location of BS is fixed like last scene, and we also additionally add regular TUAVs to densify both UL and DL. As we mentioned before, only 3 users can be connected to each TUAV in this scenario, because each user occupies 2 RBs. After that, we determine the location of TUAVs and user's association matrix as described in the previous scenario.   • Special TUAVs (Assisting DL): In this scenario, the only difference between last scene is that TUAVs can only receive messages from BS through cables, then forward them to users via wireless DL, so the TUAV can connect up to 6 users. Moreover, we also run the three algorithms mentioned above to minimize the DL exposure.
The effectiveness of different types of TUAVs in reducing EMF exposure is shown in Fig. 6 and Fig. 7. In this scenario, the number of active users K varies from 6 to 48. Among these active users, 20% of them use mobile phones mainly for making calls with required UL voice rate 5 Mbps, the remaining 80% use mobile phones for data with a high data rate requirement, i.e., 50 Mbps for UL. For DL, the rate required by the voice is 5 Mbps, while that of data is 100 Mbps. Fig. 6 reveals that using TUAV to densify DL are not an effective way to reduce EMF exposure. This is because EMF exposure from DL is almost negligible compare to that of UL. To show this more clearly, the EMF exposure separately generated by DL and UL in different scenarios are shown in Fig. 7. As can be seen, there is a great difference in order of magnitude, i.e., UL radiation is about 10 6 times that of DL. In contrast, using TUAVs that can densify UL could effectively reduce the EMF exposure, especially green TUAVs.
Because we use regular TUAV to densify both UL and DL, each user will occupy two RBs and the maximum number of users connected to TUAV will be halved. Therefore, the effect of regular TUAV is not as good as the green TUAV, given that they have the same amount of RBs which are relatively limited. Moreover, green TUAVs, with receiving only RF chains, have lower complexity than regular TUAVs.

C. Green Fixed SCs and Green TUAVs
In the following, we analyze the performance of green TUAVs and green SCs, both of which are used to densify the UL. The only difference between them is that TUAVs are much more flexible, while SCs have a fixed deployment, and they are evenly distributed in the area. Since the achieved UL data rate can be less than the required one due to the limitation in the devices' transmit power, we define the satisfied-users ratio as the ratio of the number of users who satisfy their rate constraint (14b) to that of total active users. We study the enhancement of the satisfied-users ratio when assisted by M = 4 green TUAVs or green SCs. The locations of SCs are fixed, while TUAVs have N = 36 GS to choose from. We also considerK = 240 residents distributed in the area, 25% of them are active, and each of them has the same data rate requirement, i.e., data users. Fig. 8 plots the satisfied-users ratio versus required UL data rate. As can be seen, when the required data rate is less than 50 Mbps, each network can meet the needs of almost all users. However, as the data rate continues to increase, there will be significant gaps between each scenario. For instance, considering a 100 Mbps required data rate, the fixed SCs can achieve only a gain in the satisfied-users ratio of 70% compared to BS only approach. On the other hand, the proposed TUAVs scheme boosts it with up to 400% compared to the BS only scenario.
Such a significant gain is mostly due to the reduced area that each gNB has to cover. Therefore, the power that the user should transmit to achieve a required data rate decreases significantly below the maximum allowable transmit power, escalating the satisfied-users ratio.
Furthermore, Fig. 9 depicts the UL EMF exposure versus the required UL data rate. This figure shows that as the required UL data rate increases, the total EMF exposure increases up to a certain value then saturates. This can be attributed to the fact that each user can not accommodate transmit power higher than his maximum level P max , leading to a drop in his satisfied rate, as in Fig. 8, and saturation   SCs, respectively. Combined with Fig. 8, we can see that the proposed scheme can improve the satisfied-users ratio, equivalent to the UL data rate while maintaining relatively low EMF exposure.
In Fig. 10 with required UL rate of 5 Mbps. In the case of BS only, the curve increases almost linearly from the start with the average number of active users. While the other two curves of UL EMF exposure exhibit the same basic shape, the slopes are small at first, then increase to certain values and remain stable due to the limited available RB. In general, using TUAVs can adapt to changes in user distribution. For example, the proposed scheme can reduce the exposure, for K = 50, by 65% and 30% compared to fixed BS only and SCs, respectively.
In Fig. 11, we consider a scenario ofK = 240 residents, and 40% of them are active users and divided into two categories as scenario of Fig. 10. After that, the solid lines represent radiation with green TUAVs assistance, and the dashed line parallel to the x-axis represents that of fixed deployed SCs. Moreover, the number of TUAVs and SCs M ∈ {2, 4, 8}, while that of GSs N spans from 2 to 128. As can be seen from Fig. 11, increasing the number of TUAVs or SCs is a useful way to reduce EMF exposure, but it's also costly. With the same number of TUAVs or SCs, TUAVs provide greater performance than SCs, even with M = N , i.e., also a fixed deployment for  also N = 36 evenly distributed GSs andK = 240 residents throughout the region, 25% of them are active, i.e., K = 60 users. We consider the exposure metric in terms of whole body SAR with two levels: i) the ICNIRP limit, i.e., SAR limit = 0.08 W/kg; ii) a stricter limit, i.e., SAR limit = 0.0016 W/kg, with a further safety factor of 50 compared to the ICNIRP limit. In Fig. 12, we show the average UL rate when the whole body SAR limit spans from 2 × 10 −8 W/kg to 2 × 10 −1 W/kg. Moreover, we change the value of α n , from 2 to 3.5, while keeping α l = 2, to see the impact of the path loss exponent for NLoS on the performance. As shown in Fig. 12, the average UL rate would go down if the path loss exponent α n goes up, but the change of α n does not affect the trend of the curves. Besides, it can be seen from the figure that the average rate increases as the threshold are relaxed. However, when the limit is increased to around 2 × 10 −3 W/kg, further relaxation of the threshold does not improve the average rate. This can be explained by the hardware limitation of the user's transmit power. Moreover, the use of TUAV improve the average rate by about 50 Mbps compared to fixed SCs at ICNIRP and stricter SAR limits.
In Fig. 13, the maximum transmit power changes from 20 dBm to 36 dBm, while the SAR limit remains the same at ICNIRP limit 0.08 W/kg or strict limit 0.0016 W/kg. As can be seen that the addition of auxiliary green SCs or TUAVs can effectively increase the average rate. Besides, when the ICNIRP limit is considered, the average rate is positively correlated with the logarithm of the power threshold. When the SAR limit is relatively strict, it can be seen that the curve saturates after rising, which can be explained by (21), i.e., the power allocation policy is limited by the exposure rather than the hardware constraint in this region. Moreover, under the same case, the TUAVs can adapt to the distribution of users, thus improving by an additional 50 Mbps compared to fixed SCs.

VIII. CONCLUSIONS
A prevailing theory in the non-scientific community claims that EMF exposure is uncontrolled and exponentially increasing due to gNB densification. Also, there is a debate about the adverse health impacts due to the long-term non-thermal exposure to RF radiations from BSs. In order to reduce the population exposure to EMF, we proposed a novel architecture with TUAVs carrying green antennas to assist the communication system. According to relevant work and confirmed by our numerical results, the UL exposure is the dominant factor, e.g., 6 orders of magnitudes higher than the DL exposure for the considered simulation settings. In particular, we formulated two optimization problems either to minimize the UL EMF exposure and ensure a target data rate or to maximize the sum rate with an EMF constraint. Moreover, we proposed several low-complexity algorithms to solve the optimization problems. The simulation results illustrated that using green TUAV to densify UL can provide good results when the number of RBs is limited. On the other hand, comparing fixed green SC to mobile green TUAV, the latter tend to provide a better performance, both EMF exposure and satisfied-users ratio, reducing by 23% and boosting by 400% respectively when required UL data rate is 100 Mbps. The proposed scheme for cellular system design with EMF constraint can significantly improve the average UL rate by almost 15% and 350% compared to SCs and BS only architectures when considering the ICNIRP exposure limit. Therefore, we believe that the proposed architecture can be of interest for both the general public, cellular operators, and the research community, as reducing the exposure can alleviate the adverse impacts due to the EMF exposure.