Topology optimization of microwave frequency dividing multiplexers

We use material-distribution-based topology optimization to design a three-port frequency dividing multiplexer at microwave frequencies. That is, by placing a good electric conductor inside the design domain, we aim to design a passive device that splits the incoming signal’s frequencies into two frequency bands and transmits them to their respective output ports. The Helmholtz equation models the time-harmonic wave propagation problem. We solve the governing equation using the finite element method. The adjoint variable method provides the required gradients, and we solve the topology optimization problem using Svanberg’s MMA algorithm. In this study, we present a technique for modeling the distribution of a good electric conductor within the design domain. In addition, we derive a power balance expression, which aids in formulating a series of three objective functions. In each successive objective function, we add more information and evaluate its impact on the results. The results show that by selecting a suitable objective function, we achieve more than 93.7%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$93.7\,\%$$\end{document} transmission for both the frequency bands. Moreover, the numerical experiments suggest that the optimization problem is self penalized and is sensitive to the initial design.


Introduction
The advent of wireless networks, which was fuelled by breakthroughs in microwave engineering, has revolutionized telecommunication. In the last two decades, the use of wireless and connected home devices, such as cell phones, smart speakers, security cameras, and laptops has increased significantly, resulting in an exponential increase in wireless data traffic. In 2020, global data consumption increased by 30.4 % over the previous year, and it is estimated that there will be nearly 30 billion connected devices by 2023 (Cisco 2020; PwC 2021). To accommodate this increase in data traffic, future communication systems and wireless devices need more efficient and improved designs.
Frequency dividing multiplexing (FDM) plays a key role in data transmission for telecommunication (Sakai 1986;Li and Stüber 2006;Snyder et al 2015;Ma et al 2017). It enables a high level of system miniaturization by integrating sub-systems that operate at different frequency bands. On one hand, a frequency dividing multiplexer allows multiple signals of different frequencies to share a single communication channel, such as a waveguide or an antenna. On the other hand, the device divides the frequencies of the communication channel into multiple non-overlapping frequency bands and delivers them to their corresponding output ports. To prevent signal overlapping, sufficient isolation between neighbouring frequency bands is essential. In addition, a sharp roll-off between the adjacent frequency bands is needed to avoid wasting parts of the frequency spectrum. AM and FM broadcasting methods are the most common use of FDM at RF frequencies (Chaparro 2015), and communication satellites (Kharchenko et al 2015), modern radars (Sebt et al 2009), and wireless LAN (Keller and Hanzo 1996) are just a few examples of the application of FDM at microwave frequencies.
Responsible Editor: Palaniappan Ramu In the literature, various concepts are used to design FDM devices (Snyder et al 2015;Cameron and Yu 2007). Two commonly used design concepts of FDM are the cavitybased and the directional-filters approaches. The cavitybased approach relies on the modal field distribution inside a given cavity. The positions of the FDM output ports are properly selected at the boundaries of the cavity to enable efficient signals' coupling (Yassini and Yu 2015). However, the field distribution inside a loaded cavity changes based on the loading condition, and therefore, the device needs to be redesigned to account for the new loads. On the other hand, FDM based on directional filters rely on loading a transmission line with a series of resonant filters. Each filter stage is designed separately and aims to couple a specified frequency band to an output port (Cameron and Yu 2007). Employing techniques such as adding transmission zeros can improve the transition slope between the neighbouring stages (Sorocki et al 2018). However, the final device becomes bulky, and the overall insertion losses of the device increase significantly.
Topology optimization is an inverse design technique (Bendsøe and Sigmund 2004), which was originally introduced to design linear elastic structures (Bendsøe and Kikuchi 1988). It has since been successfully applied to other disciplines including fluids (Borrvall and Petersson 2002;Gersborg-Hansen et al 2005), acoustics (Wadbro and Berggren 2006;Dühring et al 2008;Lee and Kim 2009;Wadbro 2014;Bokhari et al 2021), and electromagnetics (Aage et al 2010;Hassan et al 2014b, a;Frellsen et al 2016;Wang et al 2017;Hassan et al 2018Hassan et al , 2020. Material-distributionbased topology optimization is the most commonly used approach, where a specified design domain is divided into small elements, and a design variable is assigned to each element to determine presence or absence of material. Typically, gradient-based optimization methods are employed to solve topology optimization problems, where the gradient information can be efficiently evaluated by adjoint-field methods (Giles and Pierce 2000). In this work, we adopt the cavity-based approach to design the FDM. However, here we fix the positions of the output ports and employ the materialdistribution-based topology optimization to design a threeport FDM that operates in the X-band (9.0-10.2 GHz). The multiplexing effect is achieved by distributing a metal inside the design domain which splits the input signal frequencies from port-1 into two frequency bands and transmits them to their respective output port (either port-2 or port-3). The design objective is to maximize the transmission of each frequency to its intended output port, minimize reflection to port-1, and limit signal cross-coupling between the adjacent output ports. Previous contributions on the FDMs problems discussed the use of dielectric materials as design materials (Nomura et al 2007;Frellsen et al 2016). In such studies, the interpolation between the design material (dielectric) and the background (void) is also a lossless dielectric material, which entails that the power imposed at the input port would either be reflected back to the same port or transmitted to the output ports.
In the current work, we present a method to design efficient three-port metallic frequency-dividing multiplexers. Here, the design material is a good conductor (copper), which can be treated as a lossless material at microwave frequencies. However, when interpolating between the good conductor (copper) and the void during the optimization procedure, the intermediate materials exhibit ohmic losses. Previous studies on topology optimization of lossy material dealt with one or two port devices, such as antennas and waveguide transitions (Hassan et al 2014b(Hassan et al , 2018. For multiport design problems, ohmic losses make the design problem sensitive to the selection of the objective function. The current work features a unique setup to study two main challenges: 1) handling the ohmic losses and 2) selecting a suitable objective function for a multiport design problem. We resolve the first issue by using a non-linear filtering approach coupled with a continuation scheme. In addition, we present a technique to model the material distribution of metal as a good electric conductor against a dielectric background. More precisely, we employ numerical investigations to select suitable values of the electric conductivity for the interpolation function. For the second issue, we devised the power balance of the setup, which we use to formulate a sequence of objective functions. Each objective function includes more information than the preceding one. The results demonstrate that providing useful additional information to the optimizer improves the transmission of the signal to its intended output port.

Problem description
Consider the two-dimensional three-port setup illustrated in Fig. 1 (a,b). Port-1 is the waveguide on the left, and port-2 and port-3 are the waveguides on the right. Here, D is our design domain. We denote waves travelling towards D as incoming waves and waves travelling away from D as outgoing waves. In addition, D and the three artificially truncated waveguides constitute the computational domain depicted in Fig. 1 (c). The boundary PEC represents a prefect electric conductor, and 1 , 2 , and 3 are the boundaries of the truncated waveguides. We aim to place metallic material inside D to make this setup serve as a FDM. In other words, we want to optimize this setup in such a way that: (a) For frequency band-1, incoming waves from port-1 only travel to port-2, as illustrated in Fig. 1a. (b) For frequency band-2, incoming waves from port-1 only travel to port-3, as illustrated in Fig. 1b).
Page 3 of 16 106 The dimensions of this setup are l w = 50 mm , l D = 100 mm , h D = 100 mm , h w = 20 mm , and h s = 30 mm , as illustrated in Fig. 1c. The electric field E is governed by Maxwell's equation where is the permeability, is the electrical conductivity, and is the permittivity. The permittivity is related to relative permittivity r by r = ∕ 0 , where 0 is the permittivity in free space. Further, we assume that ≡ 0 , the permeability in free space. We search for time harmonic solutions by using the ansatz E(x, t) = ℜ{E(x)e i t } , where the angular frequency = 2 f , f is the frequency, i is the imaginary unit, and t is the time. This assumption yields that where k = ∕c is the wavenumber, and c = ( 0 0 ) −1∕2 denotes the speed of light. Inside waveguides, the solution of Maxwell's equations can be decomposed into the transverse electric (TE) and the transverse magnetic (TM) polarizations. Moreover, waveguides are typically dimensioned and used at frequencies where only one mode (namely the dominant mode) is supported at their ports (Pozar 2011). For such cases, it could be numerically more efficient to analyze such waveguide problems by employing 2D symmetry.
We assume that the electric field is polarized normal to the plane. That is, E = (0, 0, u) , where u refers to the electric field normal to the plane as demonstrated in Fig. 1. Also, we choose the frequency band of interest so that our waveguidebased FDM only supports the dominant TE 10 mode (Pozar 2011). Under these assumptions, Eq. (2) can be reduced to the Helmholtz equation, We assume that the metallic boundaries of device are perfect electric conductors. The boundary condition for the perfect electric conductor is Based on the 20mm waveguide width, the waveguide ports only support the first propagating mode between 7.50 GHz and 15.0 GHz . In this study, we choose frequencies such that all modes except the first are evanescent. The general solution in waveguide m is thus a combination of incoming and outgoing first-mode waves, which can be written as where A m and B m are the complex amplitudes of incoming and outgoing wave, respectively, at port-m, K = √ k 2 − ( ∕h w ) 2 is the reduced wavenumber, and (4) u = 0 on PEC .
The electric field is polarized normal to the plane. (a) Incoming waves from port-1 should travel to port-2 in frequency band-1 (b) Incoming waves from port-1 should travel to port-3 in frequency band-2. (c) Computational domain for finite element discretization.
in which the scaling factor S m is selected such that We use a first-order absorbing boundary condition (Engquist and Majda 1977) to truncate the domain at m . That is, The variational form of governing equation (3) with boundary conditions (4) and (8) reads: where U = {v ∈ H 1 ( ) | v = 0on PEC } . Having solved the above problem, we compute the complex amplitude of the outgoing waves by using

The power balance and the scattering parameters
We obtain a power balance expression by choosing v =ū (complex conjugate of u) in Eq. (9), which yields Taking the imaginary part of the above equation and dividing by K, we have The local coordinate systems, as illustrated in Fig. 1 (d), are selected such that x m = 0 on m ∀ m . Thus, by using expression (5), we find (6) g m = S m sin y m ∕h w , (7) ∫ m g m g m = 1.
Find u ∈ U, such that Inserting the above expression into Eq. (12) and using that ∫ m g m g m = 1 yields where the last term represents the normalized ohmic losses. According to power balance (14), the power of the incoming waves is equal to the sum of the power of outgoing waves plus the ohmic losses.
Furthermore, it is common to use scattering parameters, S mn , to define reflection, cross-coupling, and transmission, where m denotes the port of outgoing waves and n denotes the port of incoming waves. The following expression relates the incoming and outgoing wave amplitude to the scattering parameters From power balance (14), we know that the power of an incoming and an outgoing wave is proportional to the square of complex amplitude. That is, we are interested in extremizing P mn = | | S mn | | 2 . The following expression computes P mn under the assumption that the device is only excited at port-n and that the other ports are matched We want to place a good electric conductor within D , so that the three-port setup works as two independent two-port networks. In frequency band-1, we want port-1 to only communicate to port-2, such that More precisely, if the incoming waves from port-1 belong to frequency band-1, they are transmitted to port-2, while the incoming waves from port-2 are transmitted to port-1 due to reciprocity. Furthermore, in frequency band-2, we want port-1 to only communicate to port-3, such that More precisely, if the incoming waves from port-1 belong to frequency band-2, they are transmitted to port-3, while the incoming waves from port-3 are transmitted to port-1 due to reciprocity. Henceforth, we refer to P as the power matrix.

Finite element discretization
The finite element method is used to discretize the computational domain into a uniform mesh of square elements. We use basis functions that are continuous and bi-quadratic on each element. These basis functions are denoted by 1 , 2 , … , N , where N is the number of nodes in the finite element approximation. The complex amplitude function u and the test function v are approximated by Further, we approximate the conductivity by an elementwise constant function h . We define the N D × 1 vector

Modeling of a good electric conductor
We set the incoming amplitudes A 1 = 1 and A 2 = A 3 = 0 . Using power balance (14) and the definition (16), we have In the above expression, we have the power output at all the ports and the ohmic losses. We aim to design a lossless device via material-distribution-based topology optimization by placing either a metallic material or free space in each element in D to achieve the multiplexing effect. To obtain a well-performing multiplexer, the device should possess negligible ohmic losses. The perfect electric conductor has = ∞ and contributes zero ohmic losses. Using a very large value of , we can model a perfect electric conductor and approximate a zero Dirichlet boundary condition on the metallic boundary inside D . However, choosing a very high conductivity value is neither beneficial nor practical. Good electrical conductors, such as silver and copper have finite conductivity of = 6.3 × 10 7 S∕m and = 5.96 × 10 7 S∕m , respectively, at microwave frequencies. For interpolation of material in material-distribution-based topology optimization, we need to select values of conductivity for free space ( min ) and good electric conductor ( max ) to model the material distribution inside D . To find suitable values of max and min , we conducted numerical tests with a disk and triangle-shaped electric conductor inside D , as shown in Fig. 2.
The reason for using different shapes is to experiment with structures that might contribute different amount of losses. The disk and the triangle are simple geometries that possess smooth boundaries and sharp corners, respectively. For both geometries, we vary the conductivity between two extreme values. The graphs in Fig. 2 show that for intermediate values of the conductivity, both geometries exhibit losses and that the amount of losses vary with frequency. That is, the ohmic losses are both frequency and geometry dependent. However, for both geometries, the ohmic losses are negligible for ≤ 10 −4 and ≥ 10 5 . This suggests that all the values of ≤ 10 −4 and ≥ 10 5 can model free space and a good electric conductor, respectively, with negligible ohmic losses. Hence, we choose to simulate a good electric conductor using max = 10 5 and free space using min = 10 −4 . Thus, for ≤ 10 −4 and ≥ 10 5 , We tightly choose this range for since values outside of it make the problem insensitive during the optimization. Furthermore, numerical tests reveal that the values (24) P 11 + P 21 + P 31 + Ohmic losses = 1.
outside of this range have a negligible impact on the device performance.

Design definition and filtering
We employ a non-linear filtering approach to relax the problem's self-penalization issue and impose size control on the design. More precisely, we apply the fW-mean-based filtering framework (Wadbro and Hägg 2015), and harmonic mean-based filters (Svanberg and Svärd 2013) to approximate morphological operators (erode and dilate). Here, we use an extended domain, F = D ∪ M ∪ W , illustrated in Fig. 3, to filter the design variables where M is occupied by material and W is occupied by free space. The material distribution inside F is defined by N F × 1 vector in which elements are sorted such that the elements in D come first, the elements in M comes second, and the elements in W come last.
In the discrete setting, the harmonic erode and dilate operators act on and are defined as follows: and where r is the filter radius, > 0 is a parameter that controls the non-linearity, f and f D have entries and respectively, and N F × N F weight matrix W r implicitly defines neighborhood of element i. More precisely, W r = D −1 G r . In N F × N F matrix G r , the entry g ij = 1 if the distance between x i (centroid of element i) and x j (centroid of element j) is less than or equal to r, else g ij = 0 . In addition, the diagonal matrix D = diag(G r 1 N F ) , in which all the entries of N F × 1 vector 1 N F are equal to 1.
Note that the material distribution inside M and W remains fixed during optimization, while only the material distribution inside D changes. We define a N D × 1 design vector for material distribution in D . Moreover, we form a harmonic close operator by combining the erode and dilate operators in a series, providing us a size control of radius r on free space. The harmonic close operator, which acts on , is defined as of all ones, and C r, ( ) provides a filtered version of the design vector . The element values of conductivity have components where C r, ( ) E = 1 corresponds to E = 10 −4 (free space) and C r, ( ) E = 0 corresponds to E = 10 5 (good electric conductor).

Objective function
The aim of this study is to design a frequency dividing multiplexer with properties as discussed in Sect. 3. The essential requirement is that the three-port setup should work as two independent two-port networks. That is, in frequency band-1, we want maximal transmission between port-1 and port-2, and in frequency band-2, we want maximal transmission between port-1 and port-3. To this end, we investigate three objective function formulations. Each consecutive objective function includes more information obtained by adding additional terms from power balance expression (24).
In addition, to writing the objective functions in a form inspired by the power balance expression, we also write them in the exact form that we use in the numerical experiments. In the numerical experiments, we use the method of moving asymptotes (MMA) by Svanberg (1987) to solve the optimization problem. MMA provides multiple alternatives for solving an optimization problem, including standard minimization and least square problems. The reason that we provide the exact formulation of the optimization problem that enters the optimizer is that this may significantly affect the performance-both in terms of number of iterations required to find a minimizer, but also in terms of the performance for the optimized design. Henceforth, we denote the formulation we use in MMA as the "MMA formulation" of the optimization problem. (30) Finally, before stating our three objective functions, we remark that all of these objectives push the power matrix for the optimized design towards form (17) in frequency band-1 and form (18) in frequency band-2. This is further expanded in Remark 2 below.

First optimization problem
For the first objective function, we only optimize with information on transmission. More precisely, the objective is to maximize P 21 for frequency band-1 and P 31 for frequency band-2. Due to reciprocity, we also implicitly maximize P 12 for frequency band-1 and P 13 for frequency band-2. For the first study, we formulate the following objective function to be minimized: where Q 1 and Q 2 denote the number of frequencies in frequency band-1 and frequency band-2, respectively.

Remark 1
The reason that we use 1 − P 21 and 1 − P 31 is that the ideal (maximum possible) value for these normalized powers is 1. So each term in the sums above is positive, which paves the way for a straightforward transition to the MMA formulation.
In each iteration of topology optimization, we take the following steps to determine power P m1 for each wavenumber 1. For filtered design vector C r, ( ) , wavenumber k q , and incoming wave amplitude A 1 , solve finite element problem (19) to compute u. 2. Compute the amplitude of outgoing waves B m using Eq. (22). 3. Use expression (16) to compute power P m1 , m = 1, 2, 3.
Here, we aim to minimize J 1 and hence solve the following optimization problem Remark 2 An ideal design for the FDM with objective J 1 has P 21 = 1 in frequency band-1 and P 31 = 1 in frequency band-2. For frequency band-1, P 21 = 1 together with the power balance implies that P 11 = P 31 = 0 , min ∈A J 1 (C r, ( )), this together with reciprocity implies that P 13 = 0 and P 12 = 1 , which (again by a power balance law) implies that P 22 = P 23 = 0 , which by reciprocity implies that P 32 = 0 . An analogous argument yields that P 31 = P 13 = 1 and P 11 = P 21 = P 12 = P 32 = P 33 = 0 in frequency band-2. Thus, an ideal design that perfectly minimizes J 1 has a power matrix of form (17) in frequency band-1 and of form (18) in frequency band-2.
The standard minimization formulation in MMA is suitable for solving problem (33) where we minimize the difference to maximize P 21 and P 31 . Hence, we do not use the quadratic term provided in the MMA formulation. To solve our optimization problem, we use the following MMA formulation where Q = Q 1 + Q 2 denotes total number of frequencies, and the Q × 1 vector y represent artificial optimization variable that aids in the formulation of optimization problem. We provide a separate function h i to the MMA algorithm for each wavenumber.

Second optimization problem
Only information on transmission is provided to the optimizer in J 1 . A possible aid to the objective function would be to provide information on cross-coupling or reflection, which is not present in J 1 . After optimizing with J 1 , we noticed that the values of cross-coupling are higher than that of reflection. Therefore, we include information on cross-coupling in the second objective function. More precisely, in addition to maximizing transmission, we minimize P 31 for frequency band-1 and P 21 for frequency band-2. This allows us to control another parameter while also providing the optimizer with additional sensitivities to minimize the cross-coupling. To do so, we formulate the following objective function Similarly, due to reciprocity, we implicitly minimize P 13 for frequency band-1 as well as P 12 for frequency band-2. Hence, we solve the following optimization problem To solve optimization problem (37), we obtain the MMA formulation shown below where the vector y is 2Q × 1 , which is twice the size of problem formulation (34). Thus, for each wavenumber, we provide the MMA algorithm two h i , one for maximizing transmission and the other for minimizing cross-coupling. More precisely, we provide h i as defined in Eq. (35) for i = 1, … , Q and

Third optimization problem
According to Eq. (24), we can also add information about reflection to the objective function, which was not included in the second objective function. Therefore, in the third objective, we also provide information on reflection to the optimizer. More precisely, in addition to maximizing transmission and minimizing cross-coupling, we also minimize P 11 for both the frequency bands. To do so, we define the following objective function Hence, we solve the following optimization problem To solve the above optimization problem, we obtain the following MMA formulation (37) min ∈A J 1 (C r, ( )) + J 2 (C r, ( )).
where the vector y is 3Q × 1 , which is three-times the size of the problem formulation (34). To solve optimization problem (42), we provide h i as in Eq. (35) and Eq. (39) for i = 1, … , 2Q , and In this case, we provide the MMA algorithm three h i for each wavenumber, first for maximizing transmission, second and third for minimizing cross-coupling and reflection, respectively.

On some relations between the objective functions
Also for the second and third objective, one can make the same argument as that in Remark 2. So, the ideal design for any of the objective functions in this section has a power matrix of form (17) in frequency band-1 and of form (18) in frequency band-2. The key in this argument for all objective functions is that P 21 = 1 in frequency band-1 and P 31 = 1 in frequency band-2. On the other hand, having P 11 = 0 = P 31 = 0 in frequency band-1 and P 11 = 0 = P 21 = 0 in frequency band-2 does not imply that the design is ideal-it is, for example, possible that all the incoming energy gets burned due to the ohmic losses.
Another, possibly interesting relation between the partial objectives can be found by rearranging power balance (24). It holds that Thus, for lossless design the third objective function is a factor 2 times the first objective function. However, for design that exhibit losses, the partial objectives are not perfectly aligned.

Sensitivity analysis
This study aims to maximize or minimize P mn , which depends upon the complex amplitudes, A m and B m , of the incoming and the outgoing waves, respectively, at m . The values of A m are fixed, therefore we can only influence B m to optimize P mn . We employ the adjoint-based method for computing the gradient of B m with respect to electrical conductivity, E .
Let be a perturbation of the conductivity. The corresponding first-order perturbation of the outgoing complex amplitude B m (22) is where u is the first-order perturbation of u.
Similarly, since the system matrix A is linear in , the perturbed Eq. where M E is a mass matrix of element E in D and has entries u E and (z m ) E are sub-vectors of u and z m , respectively, which correspond to nodal values in element E .
By differentiating Eq. (31) and using the chain rule, we obtain the gradient of B m with respect to design variable E , as presented by Hägg and Wadbro Hägg and Wadbro (2017).

Continuation approach
Previous studies by Hassan et al (2014bHassan et al ( , 2015 have shown that the problem of designing a metallic antenna to maximize transmission is self-penalized. Based on power matrices of form (17) and (18), we choose to maximize transmissions P 21 for frequency band-1 and P 31 for frequency band-2, respectively. Maximizing the transmission implicitly minimizes the ohmic losses, which makes the optimization problem self penalized. Thus, no explicit penalty is required to suppress the intermediate values of design variables. As stated earlier, to solve this issue, we employ a continuation approach for the filtering of the design variables. That is, we solve a sequence of optimization problems, each starting with the previously optimized solution as initial guess. At each continuation step, we decrease the non-linearity parameter when the residual norm of the KKT (or the first-order optimality) condition is smaller than 10 −3 . More specifically, we use = 10 2− n 2 where n = 0, 1, … , 12 . As we decrease , the filter's non-linearity increases. The blurring caused by the filtering implicitly imposes losses. The non-linearity of the filter increases with each successive step, which is associated with and negatively correlated to the amount of induced losses. Thus, in this case, the non-linearity of filter reduces blurring in the design in addition to providing size control and mesh-independent design.

Numerical results
We use a set of 11 frequencies in each band for the topology optimization of the three port setup. More precisely, in frequency band-1, we use the frequencies 9.0, 9.02, … , 9.2 GHz , while in frequency band-2, we use the frequencies 10.0, 10.02, … , 10.2 GHz . Using the definition of wavenumber k and angular frequency w, we compute k q for each frequency in both bands. The design domain D is discretized into 320 × 320 elements. The wave propagation here is well resolved, with ∕ x ≈ 106 elements per wavelength at the lowest frequency 9.0 GHz and ∕ x ≈ 94 elements per wavelength at the highest frequency 10.2 GHz , where is the wavelength and x is the element size of finite element discretization. The mesh resolution in this case is determined by the desired design resolution in D and not by the wave propagation.
For all experiments, the final optimized designs are evaluated using a real electric conductor copper with = 5.96 × 10 7 S∕m over the frequency range 8.8-10.4GHz with an increment of 0.01 GHz . Recall that the first objective function measures the transmission for both the frequency bands. Therefore, we evaluate J 1 to compare all the optimized designs. As a reference value, J 1 = 12.2225 for the free space (empty cavity) initial design. That is, C r, ( ) E = 1 , ∀ E in D .

Maximizing transmission
For the first study, we employ our basic objective function J 1 formulated in Eq. (32). This objective function maximizes the transmission P 21 for frequency band-1 and P 31 for frequency band-2. The optimized design and frequency response are shown in Fig. 4. Within frequency band-1, more than 93.5% of the power is transmitted to port-2, with an average transmission of more than 95.1% . In addition, the reflection and the cross-coupling are less than 3.2% and 4.4% , respectively, for all frequencies in the range 9.0 − 9.2GHz . Similarly, within frequency band-2, more than 92.5% of the power is transmitted to port-3, with the average transmission to port-3 is more than 93.9% . In addition, the reflection and the cross-coupling are less than 0.7% and 7% , respectively, for all frequencies in the range 10.0 − 10.2GHz . Furthermore, J 1 = 1.1987 at the end of optimization implies that the transmission can further be improved by minimizing the cross-coupling and the reflection terms. At the completion of optimization, the element entries of filtered design vector C r, ( ) are nearly binary. More precisely, we employ the measure of non-discreteness ( M nd ) proposed by Sigmund (2007), which is For this measure, M nd = 0 % means that the element values of the filtered design vector are binary (either 0 or 1), whereas M nd = 100 % means that the design is totally gray and all the element values are 0.5. M nd = 0.0213 % for the design presented in Fig. 4, which indicates that the optimized design is essentially lossless. Power balance states that sum of power is always equal to 1. Here, we verify that the design is essentially lossless by summing the outgoing powers P 11 + P 21 + P 31 . This holds for all the numerical results presented in this work.

Maximizing transmission and minimizing cross-coupling
As stated previously, the coupling values are higher than that of reflection in the first study. Therefore, in the second study, we maximize transmission as well as minimize cross-coupling by employing the second objective function defined in expression (37). The optimized design along with its frequency response is shown in Fig. 5. Within frequency band-1, more than 93.7% of the power is transmitted to port-2, with an average transmission of more than 96.1% . Furthermore, the reflection and the cross-coupling is less than 1.7% and 4.1% , respectively, for all frequencies in the range 9.0 − 9.2 GHz . Similarly, within frequency band-2, more than 94.4% of the power is transmitted to port-3, with an average transmission of more than 95.2% . Furthermore, the reflection and the cross-coupling are less than 0.9% and 2.3% , respectively, for all frequencies in the range 10.0 − 10.2 GHz . As compared to the first study, the transmission is higher and the cross-coupling is lower for both the frequency bands. In addition, the optimized design exhibits a value J 1 = 0.9583 , which indicates that the second objective function is better than the first objective function, and hence the design provides a better multiplexing effect. Furthermore, as shown in the frequency response (Fig. 5), the multiplexing effect extends well beyond the target frequency bands, with more than 90% of power transmission to port-2 for frequencies in the range 8.9 − 9.4 GHz and more than 88% of power transmission to port-3 for frequencies in the range 9.9 − 10.3 GHz.

Maximizing transmission, minimizing cross-coupling, and minimizing reflection
Along with optimizing transmission and cross-coupling, the third objective function formulated in expression (41) also aims to minimize the reflection term P 11 for both frequency bands. The optimized design along with its frequency response is shown in Fig. 6. The optimized design obtained from the second and third objective function are almost identical. Similarly, the frequency responses are quantitatively similar. In addition, J 1 = 0.9524 for the third study. This is comparable to the second study, implying that adding the reflection term to the objective does not provide a significant benefit in this case. However, previous studies have shown that including the reflection 8.8 9 9.2 9.4 9.6 9.8 10 10

Validation of results
We use COMSOL Multiphysics to validate the results of our numerical experiments. For the validation, we select the final design from the second objective function depicted in Fig. 5. A MATLAB code that generates the coordinates was used to trace the boundaries of optimized design. Using these coordinates, the optimized design was created in COMSOL Multiphysics. We use triangular elements with a maximum side length of 1 mm for the discretization in COMSOL Multiphysics. The frequency response computed by our algorithm in MAT-LAB is compared to the frequency response computed in COMSOL Multiphysics. The two computations exhibit a perfect match over the simulated frequency band and thus validate our results, as shown in Fig. 7(c). In addition, we show in Fig. 7(a,b) the electric field distribution Re(u) at 9.0 and 10.0GHz to demonstrate the multiplexing effect of the optimized device.

Sensitivity to initial design
For the designs presented in Figs. 4, 5, and 6, we initiate the optimization process with free space inside D , and the optimizer only places material close to the boundary walls, with no material in the center of D . Therefore, it is reasonable to investigate the optimization problem with different initial designs. To study this further, we perform two more numerical experiments with different initial distribution of material inside D . In the first experiment, we use a circleshaped initial design, and in the second, we use a triangleshaped, as illustrated in Fig. 8. We use both the second and third objective functions in these numerical experiments. However, the value of J 1 is comparable for both the objective functions, which indicates that for this problem setup, adding information on reflection does not significantly improve the multiplexing effect. Therefore, in this section, we only present the results with the second objective function. Electric field distribution Re(u) is plotted: (a) at 9.0 GHz where the incoming waves from port-1 travel to port-2 and (b) at 10 GHz where the incoming waves from port-1 travels to port-3.
Design validation: (c) Frequency response computed by our algorithm implemented using MATLAB compared to frequency response computed using COMSOL Multiphysics Figure 8 shows the optimization results using free space, circle-shaped and triangle-shaped initial designs. For the case with free space initial design, the optimization algorithm initially places gray material close to the boundaries of D but not in the middle. Moreover, the final design we acquire has no material in the middle of D . For cases with circle-shaped and triangle-shaped initial designs, we obtain gray material close to the boundaries and the center of D . The material in the middle changes shape but does not disappear during the optimization process in both cases. This investigation highlights how the problem is affected by the initial design. Within frequency band-1, we obtain more than 90% transmission to port-2 for all three initial designs, with the circle-shaped initial design providing a flat response and nearly 100% transmission to port-2. Moreover, the difference between J 1 = 0.5370 for circle-shaped initial design and J 1 = 1.5570 for triangle-shaped initial design suggests that the former provides a much better multiplexing effect.  However, for some frequencies in frequency band-2, the transmission to port-3 falls below 90% for circle-shaped and triangle-shaped initial designs.

Discussion
In section 4.4, we present three MMA formulations for our optimization problems, and explain how we increase the number of functions h i for each objective function. We can, however, decrease the number of functions provided to the algorithm, which affects the design and its performance. We decreased the number of functions by using the same problem formulation as in (34), where instead of providing one h i for each term to the algorithm, we just provide two, one for each frequency band. More precisely, for the first objective function, we use the following MMA formulation Similarly, we provide only two h i for the second and third objective function as well. With only two h i provided to the MMA, J 1 = 1.3952 , 1.0545, and 1.0235 for the first, second, and third objective function, respectively, with free space initial design. Providing information on cross-coupling and reflection also improves multiplexing in this case. Moreover, with the second objective function, J 1 = 7.9185 and 3.0855 for the circle and triangle-shaped initial designs, respectively, by supplying only two h i to the MMA. These designs are inferior compared to the designs presented in section 5. In short, we want to highlight two key points from these numerical experiments. First, the results suggest that in this case, providing more information (cross coupling and reflection terms, as well as their gradients) aids the optimizer in finding a better solution, resulting in an improved multiplexing effect. Second, providing information to the optimizer as separate functions-as is the standard practice for least squares problems, where specialized algorithms, such as Gauss-Newton and Levenberg-Marquardt, make use of this information-yields better performing optimized designs. In all three objective functions, we use an explicit term to maximize the transmission because the minimization (53) min ,y of the reflection and cross-coupling terms does not imply the maximization of the transmission term. This is due to ohmic losses, which can be seen in Eq. (24). According to the system power balance, a possible solution to minimize the reflection and the cross-coupling is to increase ohmic losses. As a result, in a wave propagation problem with losses, an explicit term is required to maximize transmission. Maximizing transmission makes the optimization problem self-penalized and the final designs are essentially black and white without using an explicit penalty. More precisely, the maximum measure of non-discreteness ( M nd ) for all optimized designs in this paper is 0.0223 % . In a lossless case, however, we can implicitly maximize transmission by minimizing cross-coupling and reflection (Wadbro 2014;Bokhari et al 2021).
We experimented with numerous initial designs in this study, including square-shaped initial designs, various sizes of circles, and triangular forms. We also experimented with different values of gray material initial designs. The results show that many good designs are possible with various initial designs and there are multiple good local minima. Here, we present three initial designs to optimize our three-port setup, all of them provide different designs and produce a good multiplexing effect. However, we only achieve a flat response with nearly 100% transmission in frequency band-1 for the circle-shaped initial design. This demonstrates that the initial design has an impact on the final optimized design, and we typically end up at a local minimum.

Conclusions
FDM plays a crucial role in data transmission for communication systems. Because of the complexity of the design task, the design of FDM devices is typically carried out manually by highly experienced engineers. To automate the design procedure, we present a topology optimization method relying on an accurate numerical solutions of the wave propagation. Using metals as design materials in topology optimization of electromgentic problems results in ohmic losses associated with intermediate designs. The ohmic losses and the multiport features of the FDM's setup make the optimization problem sensitive to the formulation of the objective function and the solution strategy. The sole use of the maximization of power transmission, as an objective function, reveals a self-penalization of such formulation: the maximization of the transmitted power could be easily achieved by minimizing the ohmic losses in the design domain. The filtering method relaxes the self-penalization issue during the design process by controlling the ohmic losses which enables the optimization procedure to converge to satisfying solutions. The power balance of the setup, however, reveals that other factors can provide additional useful information which, when included in the objective function formulation, can result in better solutions. In other words, the system power (or energy) balance is essential to guide the formulation of objective function. Also, in designing metallic multiport electromagnetic devices, the design problem is likely to be formulated as a multiobjective optimization problem that combines various power (or energy) measures.
Acknowledgements The computations were partially performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the High-Performance Computing Centre North (HPC2N).
Funding Open access funding provided by Umea University. This work was partially supported by the Swedish strategic research programme eSSENCE and by the Higher Education Commission (HEC) of Pakistan.

Code availability
The code is available from the corresponding author upon reasonable request.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Replication of Results
All the details necessary to reproduce the results have been defined in the paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.