Efficient Solvers for Coupling-Aware Beamforming in Continuous Aperture Arrays

1 Ef ficient Solv ers for Cou pling-A wa re Beamforming in Continuo us Aperture Arrays Geonhee Lee, Kwonyeol Park, Hyeongjun Park, Jinw oo An, and Ju n il Choi Abstract —In continuous aperture arra ys (CAP As), careful consideration of the un derlying physics is essential, among which electroma gnetic (EM) mutual coup ling plays a critical role in beamf orming perf ormance. Building on a physically consistent mutual coupling model, th e beamf orming design i s f ormulated as a function al opti mization wh ose optimality condition leads to a Fredholm integral equation. Th e incorporation of the coupling model, howeve r , substantially i n creases computational complexity , necessitating efficient an d accurate integral equation solvers. In this letter , w e propose two efficient solvers: 1) a co ordinate-transform ation-based kernel a pproximation that preser ves the operator structure while alle viating discretization demands, and 2) a direct lower -upper (LU)-based solver th at stably handles the Nystr ¨ om-discretized system. Numerical results demonstrate improv ed accuracy and reduced computational ov erhead compared to con ventional methods, with the LU-based solver emerging as an efficient and scalable solution for large- scale CAP A opti mization via offline factorization. Index T erms —Beamforming, continu ous aperture array , mu tual coupling, polarization. I . I N T R O D U C T I O N Continuou s aperture arrays (CAP As) are emerging as a promising electromagn etic ( EM) arch itecture ch aracterized by aperture - lev el continuity [1]. Such continuity relaxes th e discrete spacin g con straints of conventional array s, allowing near-continuou s man ipulation o f surface cur rent d istributions. Consequently , precise amplitude and phase con trol can be achieved across the ape rture, enabling th e array ’ s EM p roper- ties to be effecti vely exp lo ited. Pr io r works have investigated various design strategies for CAP As, includ ing calculus-o f- variations-based fo rmulations to characterize energy-efficient current distributions [2] and linear be amformin g schemes to improve tractability in practical d esign [3 ]. The weig hted min- imum mean-sq uare err or (WMMSE) framework enables sum- rate o ptimization in CAP A-enab led transceivers [4]. T ogeth er , these efforts h ighlight the need fo r careful modelin g o f EM proper ties in CAP A system design. More recen tly , an EM c o upling kernel fo r CAP As and a correspo n ding numer ical me thodolog y h a ve b een proposed , laying th e foun dation for p hysically consistent CAP A beam- forming [5]. W ithin this framew ork, tw o appro aches have been propo sed to solve the co u pling-aware integral equation for G. Lee, H. Park, J. An, and J. Choi are with the School of Electrical E ngi- neering , Ko rea Adv anced Institute of Science and T echnolog y , Daejeon 34141, South Kore a (e-mail: { peppermi nt01, m ika0303 , freddy1 , junil } @kai st.ac.kr). K. Park is with the S.LSI Divisio n, Samsung Electronic s Company Ltd., Hwaseon g-si, 18448, South Kore a, and also with the School of Electrica l Engineeri ng, Korea Advance d Institute of Science and T echnol ogy , Daeje on 34141, S outh Korea (e-mail: kwon10.park@ka ist.ac.kr). optimal beamf orming. Alth ough bo th a im to solve th e equation effecti vely , each exhibits n otable limitations that affect stability or computation al effi ciency . Kernel a pproxim ation (KA) pro- vides closed- form inv erse operator s that offer physical insight and fast computatio n, wh ile explicit inv ersion c a n lead to numerical instab ility . Meanwh ile, the c o njugate gr a dient (CG) method avoids direct inv ersion throu gh an iterative proced u re. Howe ver, each iteration requir es d ense operato rs, an d the iteration c o unt g rows with a p roblem size and an o perating frequen cy , lead ing to high comp utational co st. Incorp orating mutu al coupling into existing CAP A al- gorithms fur ther intensifies th e se challen ges. For example, WMMSE-based CAP A algorithm s requ ire multiple solves o f the in tegral equation , m aking n umerical stability and ef fi- ciency essen tial. According ly , high appr oximation accura cy is required even with c o arse Gauss–Legen d re (GL) sampling . A t the same time, the need for r epeated solutions of th e integral equations necessitates low-complexity strategies capab le o f balancing accura cy , robustness, and efficiency . In this letter, we develop two efficient comp utational strate- gies for solv in g th e EM mutual coup ling p roblem in CAP As, with a fo c us o n enablin g fast an d a ccurate appro ximation while supportin g reliab le numerical beh avior . Th e main contributions of this letter ar e summarized as follows. • W e first develop a po lar-trigonometric kernel approx i- mation ( PKA) th a t enables rapid, high-fidelity mo deling while p reserving analytical structure, sup porting r e liab le coupling -aware contin uous-aper ture be a m forming . • W e also dem onstrate that a direct lower -u pper (LU)-based solver reliably resolves the Nystr ¨ om-discretize d sy stem, achieving stable perform ance with reduced compu ta tio nal complexity . • Nu merical r e sults dem onstrate that the pro posed metho ds achieve superior accu racy an d lower com putational cost than conventional KA a n d CG method s. Collectiv ely , th is work establishes practical and scalab le technique s for co u pling-aware beamfo rming, enab ling numer- ically stable and low-latency optimizatio n ev en fo r e le c tr ically large con tinuous apertures. I I . S Y S T E M M O D E L A N D P RO B L E M F O R M U L AT I O N In this section, we establish th e theore tical framework for the CAP A system. W e begin by charac ter izing the EM radiation mechanism of a co ntinuou s ape rture and the resulting mutual coupling effects. Based on these fou ndations, we formu late the ar ray gain max imization p roblem and d erive the structure of the op timal current distribution. 2 A. CAP A and Rad iated F ield W e consider a transmitter arch itecture based on CAP A, where the rad iating interface is mod eled as a con tinuous surface S . The aperture is assumed to be a rectang ular planar surface lying in the xy -plane, defined as [4 ], [5] S =  [ s x , s y , 0] T     | s x | ≤ L x 2 , | s y | ≤ L y 2  , (1) where s x and s y denote the local spatial coo rdinates. L x and L y represent the physical width and height of the apertu re surface S , resp e c ti vely . According to [6 ] , the radiated electric field E rad ( r ) ∈ C 3 at an arbitra r y o b servation poin t r is expressed as follows E rad ( r ) = Z S G ( r − s ) J ( s ) d s , (2) where G ( r − s ) denotes the free-space dyadic Green’ s function and J ( s ) r epresents th e comp lex e lec tric surface curren t den - sity at location s . As a linear oper a tor , the d yadic Green’ s function maps th e current so u rce to the re sulting radiated electric field and is d efined as G ( s ) = − j k 0 Z 0  I 3 + 1 k 2 0 ∇∇  e j k 0 k s k 2 4 π k s k 2 . (3) Here, λ is the wa veleng th, Z 0 = 120 π is the intrin sic impedanc e of fr ee space, and k 0 = 2 π / λ is the wavenumber . W e consider a n a r rowband single-carrier system with a purely y -dire c ted surface curr ent over the transmit aperture, where the sourc e curren t density is expressed as J ( s ) = x ( s ) u y with u y = [0 , 1 , 0] T denoting the unit vector along the y -axis and x ( s ) representin g th e scalar current excitation at location s . Let y ( r ) denote the effecti ve electric field at the receiver located at r . Und e r line-of -sight (Lo S) p ropaga tio n and perf ect polarization alignm e nt be tween the transmitter an d a single-an tenna receiver , y ( r ) is giv en b y y ( r ) = u T y E rad ( r ) = Z S u T y G ( r − s ) J ( s ) d s = Z S h ( r , s ) x ( s ) d s , (4) where h ( r , s ) = u T y G ( r − s ) u y . In the far -field r egime, with R 0 = k r k 2 greatly exceeds the ap erture dimensions, the spherical wav efront can b e appro ximated by a plan e wa ve [5]. Under th is appro ximation, the channel resp onse depends o nly on s fo r a fixed re c ei ve lo cation r and re d uces to h ( s ) = α ( κ r ) e − j κ T r s , (5) where κ r = k 0 [sin θ cos φ, sin θ s in φ ] T with the ele vation and a zimuth a ngles of the pr opagation dir ection θ and φ , respectively . The complex gain α ( κ r ) is given by α ( κ r ) = − j k 0 Z 0 e j k 0 R 0 4 π R 0  1 − sin 2 θ s in 2 φ  . (6) B. Mutual Coupling in CAP A T ransmitter s Mutual cou pling in CAP As refers to EM intera ctions whereby the power need ed to su p port a curre nt distribution is affected by bo th its self-gen e r ated fields and the fields from neigh boring p arts of the apertur e. Based o n the ph ysical model e stablished in [5], the total electric field at the aper ture surface E tot ( s ) is decomp osed into dissipative and radiative compon ents as E tot ( s ) = E diss ( s ) + E rad ( s ) . (7) The dissipated field, accounting fo r surface resistance Z s , is modeled as E diss ( s ) = Z s J ( s ) . Conseq uently , the total field E tot ( s ) is given by [ 5] E tot ( s ) = Z S { Z s δ ( s − z ) I 3 + G ( s − z ) } J ( z ) d z . (8) The average EM transmit power P em is then e valuated by the following qu adratic functional, expressed as [ 7] P em = 1 2 Z S Z S J H ( s ) × Re { Z s δ ( s − z ) I 3 + G ( s − z ) } J ( z ) d z d s . (9) For a purely y - directed surface current J ( s ) = x ( s ) u y , th e power fun c tio nal re duces to the scalar fo r m P em = 1 2 Z S Z S x ∗ ( s ) c ( s − z ) x ( z ) d z d s , (10) where the scalar co upling kern el c ( s ) is given by c ( s ) = Z s δ ( s ) + Re  u T y G ( s ) u y  = Z s δ ( s ) + k 0 Z 0  1 + 1 k 2 0 ∂ 2 y  sin( k 0 k s k 2 ) 4 π k s k 2 . (11) W e defin e th e seco nd term as the rad iation coupling compo- nent, such that c ( s ) = Z s δ ( s ) + c rad ( s ) . C. Pr oblem F ormulation and Optimal Beamfo rming Structur e Follo wing the f ramew ork established in [5], we optimize the continuous curren t distribution x ( s ) to maximize th e array gain at the target rece iver subject to a tota l transmit power constrain t P t . This leads to th e following function al optimization pro b lem: (P1) : max x ( s )     Z S h ( s ) x ( s ) d s     2 s . t . 1 2 Z S Z S x ∗ ( s ) c ( s − z ) x ( z ) d z d s ≤ P t . (P1-a) By applying the ca lcu lus of variations to (P1), the optimal current distribution x opt ( s ) is proportion al to an auxiliary function v ( s ) , which satisfies the Fr edholm in tegral equ ation: Z S c ( s − z ) v ( z ) d z = h ∗ ( s ) , ∀ s ∈ S . (12) The optimal distribution is then given by x opt ( s ) = s 2 P t R S h ( z ) v ( z ) d z v ( s ) , (13) where scaling ensur es th e transmit power co nstraint. Hence, the key compu tational task is to solve (1 2) to recover v ( s ) , as it directly determin es the o ptimal curr e nt distribution x opt ( s ) . Howe ver, the con tinuous fo rmulation rend ers the integral eq uation generally intracta b le, n ecessitating efficient solvers. Th e next section presents n umerical strategies for th is problem , focusing on acc u rate and co mputationa lly efficient solution methods. 3 I I I . C O U P L I N G I N T E G R A L E Q UA T I O N S O L V E R S This section b r iefly revie ws conventional algorithm s as a baseline and prop oses tw o novel approaches: a po- lar–trigon ometric quadr ature for accura te KA and an LU decomp o sition metho d fo r the Nystr ¨ om-discretized system. A. Review of the Co n ventional Kernel A ppr o ximation As described in [5], the KA method aims to obtain a tractable represen tation of the cou pling kern el c ( s ) that ena bles a closed-for m inverse. Specifically , th e rad iation mu tual c o u- pling kernel c rad ( s ) is expressed in the wav enumb e r doma in throug h the two-dim ensional ( 2 D) Four ier tr a n sform, yielding C rad ( κ ) , which is g i ven by C rad ( κ ) =    Z 0 (1 − κ 2 y /k 2 0 ) 2 √ 1 −k κ k 2 2 /k 2 0 , k κ k 2 ≤ k 0 , 0 , k κ k 2 > k 0 , (14) where th e wa ven u mber vector is define d a s κ = [ κ x , κ y , 0] T . The spatial kernel c rad ( s ) is then recovered via the inverse 2D Fourier tr a n sform as c rad ( s ) = 1 (2 π ) 2 Z k 0 − k 0 Z + √ k 2 0 − κ 2 x − √ k 2 0 − κ 2 x C rad ( κ ) e j κ T s d κ . (15) T o discretize this integral, [5] em ployed Gauss–Legen dre (GL) quadr ature in Cartesian coordin a te s, where th e sampling points are co nfined to th e disk k κ k 2 ≤ k 0 . Since th e integra- tion is two-dimension al, the appr o ximation takes the form of a dou b le summ ation: c rad ( s ) ≈ M X n =1 M X m =1 a nm e j κ T nm s . (16) T o derive the quadr ature c o efficient a nm , first let ( θ n , ω n ) represent the Gauss–Legendre no de–weigh t pa ir s on [ − 1 , 1] . T o con struct wavenumber samples ( κ x,n , κ y ,nm ) within the disk k κ k 2 ≤ k 0 , we discr e tize the x -compo nent as κ x,n = k 0 θ n , w x,n = k 0 ω n , (17) where κ x,n ∈ [ − k 0 , k 0 ] is the wavenumber sample with quadra tu re weigh t w x,n . For each κ x,n , the range of κ y ,nm is giv en by h − p k 2 0 − ( κ x,n ) 2 , p k 2 0 − ( κ x,n ) 2 i . W ith a seco nd GL quadr ature, the y-comp onent is sampled as κ y ,nm = p k 2 0 − ( κ x,n ) 2 θ m , w y ,nm = p k 2 0 − ( κ x,n ) 2 ω m . The asso- ciated quadr ature coefficient is then given by a nm = w x,n w y ,nm (2 π ) 2 C rad ( κ nm ) , (18) where κ nm = [ κ x,n κ y ,nm 0] T . By utilizing (1 1) an d (16), c ( s ) can be approx imated as c ( s ) ≈ Z s δ ( s ) + M X n =1 M X m =1 a nm e j κ T nm s . (19) In prin ciple, (1 2) can be solved via the inv erse kernel c − 1 ( s ) , defined by the ortho gonality conditio n R S c − 1 ( z ′ − s ) c ( s − z ) d s = δ ( z ′ − z ) [ 5 ]. T his yields the formal solution v ( z ) = R S c − 1 ( z − s ) h ∗ ( s ) d s . When applied to the kerne l c ( s ) , the in verse kernel admits the representation c − 1 ( z − s ) = 1 Z s δ ( z − s ) + M X n =1 M X m =1 b nm e j κ T n z − j κ T m s , (20) where b nm are chosen to satisfy the orth ogonality co ndition. Therefo re v ( s ) is giv en b y v ( z ) = 1 Z s h ∗ ( z ) + M X n =1 M X m =1 b nm e j κ T n z  Z S e − j κ T m s h ∗ ( s ) d s  . (21) This closed-form app roximation provides a tractable inverse operator for scalable b eamform in g design in CAP A systems. B. Pr oposed Method I: P olar-T rigon ometric Kernel Appr oxi- mation While the Cartesian GL quadratu re in (18 ) is straightfor- ward, it become s numerica lly unstable n ear the bo undary of the pr o pagation region. As k κ k 2 → k 0 , th e d enominato r in (14) appro aches zero, causing C rad ( κ ) to diver ge. Althoug h the case κ y = k 0 remains finite d ue to the simultane o us vanishing of num erator, the integrand is o therwise unbou nded near the bou ndary , leading to large integration er rors and requirin g prohib iti vely large M for accu rate ev alu ation. T o m itigate this numerica l instability , we adopt a polar coordin ate transfo rmation com bined with a trig onometr ic sub- stitution, expressing the wa venumber c o mpone n ts as κ x = k 0 sin θ cos φ, κ y = k 0 sin θ sin φ, (22) where θ ∈ [0 , π / 2] and φ ∈ [0 , 2 π ] . Substituting these expressions in to the deno minator yields q 1 − k κ k 2 2 /k 2 0 = p 1 − sin 2 θ = co s θ . (23) Since the Jaco bian of this transformation is k 2 0 sin θ cos θ , the cos θ ter m in the deno minator is exactly canceled, th ereby removing the singular behavior when k κ k 2 → k 0 . Un der the polar–trigonom etric tr ansformatio n, the coefficients a nm are modified as ˜ a nm , given by ˜ a nm = Z 0 k 2 0 2(2 π ) 2  1 − sin 2 θ n sin 2 φ m  sin θ n w θ ,n w φ,m . ( 24) Let { ( ξ n , ω n ) } M n =1 denote the GL no des a n d weights on the referenc e in te r val [ − 1 , 1] . For the ele vation angle, the no des θ n and weigh ts w θ ,n are obtain ed by linearly mappin g this interval onto [0 , π / 2] : θ n = π 4 ( ξ n + 1) , w θ ,n = π 4 ω n , n = 1 , . . . , M . (25) Like wise, the azimuth nod es φ m and weights w φ,m are mapped onto [0 , 2 π ] as φ m = π ( ξ m + 1) , w φ,m = π ω m , m = 1 , . . . , M . (26) By combin ing (2 2), (25) , and (26), the sam pled wa venumber vectors κ nm for the pro posed method are obtained as κ nm =  k 0 sin θ n cos φ m k 0 sin θ n sin φ m 0  T . (27) This sam pling en sures that the M 2 points fo llow the Ga uss– Legendre rule over the disk, while th e trigo nometric substitu- tion handles the kernel singu larity a t θ n = π / 2 . 4 C. Pr opo sed Method II: Dir ect Nu merical solve via LU De- composition According to (11) and ( 12), th e optimal auxiliary f unction v ( s ) is o btained by solving th e integral e quation: Z S c rad ( s − z ) v ( z ) d z + Z s v ( s ) = h ∗ ( s ) , ∀ s ∈ S . (28) The Ny str ¨ o m metho d, a well-established techniqu e for solv in g integral equations [8 ] , is em ployed to enable numerical com- putation. Specifically , we discretize the spatial d o main S u sing GL quad r ature. L e t { s n , w n } N n =1 denote the q u adrature n odes and weights, wher e N = M 2 is the total nu mber of samples obtained from an M × M g rid. Applying this discretization conv erts (2 8) in to a finite-dimen sional lin ear system: N X m =1 c rad ( s n − s m ) v ( s m ) w m + Z s v ( s n ) = h ∗ ( s n ) , (29) where n = 1 , . . . , N . Th is can be expressed in m atrix-vector form as Cv = b , where v = [ v ( s 1 ) , . . . , v ( s N )] T and b = [ h ∗ ( s 1 ) , . . . , h ∗ ( s N )] T . The matrix C ∈ C N × N is defined by [ C ] nm = c rad ( s n − s m ) w m + Z s δ nm , (30) where δ nm denotes the Kronecker delta. In prac tice , the c o u- pling matrix C tends to b ecome ill-condition ed, particular ly at high sampling d ensities and for small surface im pedance values Z s . As a result, dire c tly computing the inverse of C may sign ificantly mag nify pertu rbations in the d ata, leading to severe loss of numerical ac c u racy and amplified ro und-off errors. T o e n hance numerical robustness, we instead solve the linear system Cv = b using LU deco mposition with partial pivoting. Specifically , th e matrix is factorized as PC = LU , where P is a perm utation matrix , L is a unit lower triangular matrix, and U is an uppe r trian gular matrix. Th e overall p rocedur e is summarized as follows [9]: 1) LU fact orization: Compute th e factorizatio n PC = LU using p artial piv o ting to improve nume r ical stability . 2) Forward substitution: Solve Ly = Pb by ev aluating the entries of y from top to botto m. 3) Backwa rd substitution: Solve Uv = y by ev a lu ating th e entrie s o f v fro m bottom to top . This app roach av o ids explicit matrix inversion by solv in g the system th rough stable triangular opera tio ns, thereby im - proving numerical stability and accuracy . T o furth er en hance robustness, partial piv o ting is incorp orated d uring the factor- ization to control elem ent growth and r e duce sensitivity to round -off error s. Once the LU factors are comp uted, each solve requires only f orward and b ackward sub stitutio ns, resulting in efficient com putation with predictable complexity [9]. Meanwhile, [5] also propo sed a CG m ethod to solve (1 2) without explicitly fo rming the in verse kerne l. By r e formulatin g the pr o blem as a qu a dratic functional minimizatio n , th e CG method ach ie ves stable conv ergence via the iterativ e upd ates as v ( n +1) = v ( n ) + α n p ( n ) and p ( n +1) = r ( n +1) + ξ n p ( n ) , where v ( n ) , r ( n ) , an d p ( n ) denote the current iterate, residual, and co njugate search d irection, respe cti vely . T he par ameters α n and ξ n denote the step size and u pdate coefficient. As shown in Fig. 5 of [5 ], the req uired itera tion count for the CG method grows with the operating frequency , inc r eas- ing the computatio n al burden . Since each itera tion entails a dense matrix –vector mu ltiplication with O ( N 2 ) comp lexity , the overall comp utational cost becomes substantial. In con tr ast, the LU-based solver shifts mo st of th e com- putational effort to a one- tim e factoriz a tion of the coup ling matrix. After th is prep rocessing step, subseque nt solutions can be obtained with low overhead , making the ap proach well suited for r eal-time beamf o rming un der fixed h ardware configur ations while maintain ing stro ng numer ical reliab ility . I V . S I M U L A T I O N R E S U LT S In this section , we present th e numerical results to evaluate the perf o rmance of pr oposed solvers fo r CAP A beamfo rming. Follo wing [5] , we consider a rectang ular planar ap erture S with side lengths L x = L y = 0 . 5 m and set the su rface resistance to Z s = 0 . 0128 Ω . The chann el is modeled u nder a far -field Lo S assumption , wh ere the receiver is p o sitioned at a distance R 0 = 50 m w ith elev atio n θ = 0 an d azimuth φ = 0 . T o assess the effecti vene ss o f the pr oposed ap proach e s, we compare the following four solvers: 1) K A : The conv entional kernel app roximation method based on GL q uadratur e in Cartesian coo rdinates [5]. 2) CG : The conjug ate grad ient method presented in [5], serving as an iter ati ve baseline for so lv ing the Fredho lm equation. 3) PK A (Polar - T rigo nometric KA) [P roposed] : The pro- posed solver employing a polar coor dinate transforma- tion with trigo nometric substitution to KA metho d. 4) LU [Proposed] : Th e proposed num erical so lver of the Nystr ¨ om-discretized system using LU deco mposition. The p rimary perfor mance metr ic is th e nor malized array gain, defined as [5 ] G norm = 1 P t     Z S h ( s ) x ( s ) d s     2 , (31) which re p resents the ratio between the received signal power and the total tra n smit power P t . A. Appr ox imation P erformance and Numerical Sta bility In this su bsection, we assess how accurately the proposed solvers approxim ate the tr u e solution and examine their nu- merical stability in comp arison with co n ventional benchm arks. Fig. 1 shows the n ormalized arr ay g ain as a functio n of the number o f GL qu adrature points, d enoted b y M . Because th e KA and PKA alg orithms exhibit oscillatory behavior b etween ev en and od d ord ers, the results are sampled at intervals of two to im prove visual clarity . As the CG and LU cu rves are nearly indistinguishable over the entire ran ge of M , the CG curve is o mitted. The co nventional KA m ethod app roaches the true value but exhibits irr egular and oscillatory conv ergence. Th is behavior stems fr om the singu larity of the coupling kern el near the 5 10 20 30 40 50 60 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 Fig. 1: Normalized array ga in versu s Gauss–Legen dre quad ra- ture order M for L x = L y = 0 . 5 m . propag ation bou ndary , wh ich amplifies in tegration err ors under Cartesian co ordinates. As a result, dense samp ling is r equired to o btain a stable solu tio n, leading to increased co mputation al burden at higher operating frequ encies. In contra st, the p roposed PKA method rap id ly ap proxim a tes the true value, attaining the optimal arr ay g ain with substan- tially fewer quadratu re poin ts. For exam ple, at f = 8 GHz , the required o rder for KA is 44 wh ile tha t for PKA is 28, effec- ti vely shrinking the p roblem size f r om roug hly 2000 × 200 0 to 700 × 700 . Th is improvement validates the effectiveness of the po la r–trigon ometric transf o rmation in rem oving the kernel singularity . The LU solver exhib its highly stab le conver gence to the true arr ay gain once sufficient spatial discretizatio n is ap plied. This beh avior is in sharp co ntrast with the KA-based meth ods, which exhibit less regular a nd o scillatory conv ergence, thereb y highligh ting the substantial advantage of the direct solver in numeric a l stability . Since the LU factorization can be perfor med offline an d reused for different channel realization s, the direct solver is par ticu larly attrac tive for pr actical system s requirin g both accur acy and rea l- time operation . Overall, the resu lts con firm that the pro posed meth o ds substantially im prove nu merical stability while redu cing com- putational effort compared to the conventional b e n chmark . B. Computationa l Complexity and Runtime T able I co mpares the c omputation al perfor mance of the LU solver and the iterative CG m ethod across different f re- quencies. All results are ob ta in ed with M = 25 . The LU decomp o sition is perform ed on ce an d reused for 6 25 solves with the LU solver (via forward and b ackward substitution), while the CG solver indep endently perform s 625 line a r solves. A key observation is the consistent run time of the LU- based solve, with the total execution time remain ing nea rly constant at appro ximately 1 . 9 sec despite in creasing frequ ency . This stability stems fr om the fact that th e problem matrix is deter mined solely by the apertu re geometry and operating frequen cy , enabling th e LU factor s to be pre computed and T ABLE I: Computation al complexity com parison between LU-based and CG-based solvers. Frequenc y 2 GHz 4 GHz 6 GHz 8 GHz LU decomp (1x) [ ms ] 652 .584 625.216 703.268 669.988 LU solve (total) [ ms ] 1909.334 1904.806 1901.497 1880.552 CG s olv e (total) [ ms ] 18525.238 4085 2.152 58448 .329 56571 .415 CG iteration (avg . ) 265.9 590.0 877.0 880.0 reused. As a r esult, subsequen t solves require on ly f orward and b ackward substitutions, leading to p r edictable and low latency . In contrast, the CG metho d dem o nstrates strong sensitivity to frequ ency . I ts total runtime increases fro m 18 . 5 sec at 2 GHz to over 58 . 4 sec at 6 GHz, primarily d riv en b y the growth in the average iteration co unt fro m 265 . 9 to 877 . Since e ach CG iteratio n in volves a large-scale d ense matrix– vector m ultiplication, higher freq uencies worsen the numer ical condition ing and th e r eby increase the compu tational cost. In conclusio n, th e direct LU solver delivers sup erior speed and scalability , m a king it particularly well suited for high- frequen cy CAP A beamfor m ing scenar io s. V . C O N C L U S I O N In this letter , we addressed nu merical challeng es in continuo us-apertur e b eamform ing arising from the explicit modeling of EM mutua l cou pling. W e proposed a co ordinate- transform ed kernel appro x imation that pre ser ves the analyt- ical stru cture wh ile alleviating the need for excessi vely fin e discretization. In addition , we sh owed that direct LU solver provides a reliable and n umerically stable solutio n to the Nystr ¨ om-discretized system under Gauss–Legendre q uadra- ture. Simulation results sh owed that the polar–trigon o metric kernel app roximation achieves faster an d mor e consistent conv ergence to the op timal a r ray ga in with sign ificantly fewer quadra tu re p oints than a bench mark, whereas the direct L U solver provide s h ighly r obust num e rical stability along with predictable runtime. Mor eover , offline factoriza tio n supports low-latency up dates, m aking th e p r oposed LU solver an e ffi- cient and scalable solutio n for large- scale CAP A optimizatio n. R E F E R E N C E S [1] Y . Liu et al. , “CAP A: Continuo us-aperture arrays for revo lutioni zing 6G wireless communicati ons, ” IEEE W ir eless Commun. , vol. 32, no. 4, pp. 38–45, A ug. 2025. [2] M. Qian et al. , “Cont inuous aperture array (CAP A)-based multi-group multica st communica tions, ” IEEE T rans. Commun. , vol . 74, pp. 3787– 3801, Jan. 2026. [3] C. Ouyang et al. , “Linear recei ve beamforming for CAP A systems, ” IEEE T rans. W irel ess Commun. , vol. 25, pp. 1030–1047, Jul. 2026. [4] Z. W ang et al. , “Beamforming design for continuous apertu re array (CAP A)-based MIMO systems, ” arXiv preprint arXiv:2504.0018 1 , 2025. [5] Z. W ang et al. , “Mutual coupli ng in conti nuous aperture arrays: Physical modeling and beamforming design, ” arXiv pre print arXiv:2511.11225 , 2025. [6] A. Poon et al. , “Degre es of freedom in multiple-a ntenna channels: a s ignal space approac h, ” IEEE T rans. Inf . Theory , vol. 51, no. 2, pp. 523–536, Feb . 2005. [7] A. Pizzo et al. , “Mutual coupli ng in holographic MIMO: Physical modeling and informat ion-theor etic analysis, ” IEEE J. Sel. Are as Inf . Theory , vol. 6, pp. 111–126, May 2025. [8] R. Kress, Linear Inte gral Equations , 3rd ed. Springer , 2014. [9] G. H. Golub et al. , Matrix Computat ions , 4th ed. Johns Hopkins Uni versi ty P ress, 2013.