Ali Momeni

Babak Rahmani

Matthieu Malléjac

Philipp Del Hougne

Romain Fleury
email: romain.fleury@epfl.ch

Backpropagation-free Training of Deep Physical Neural Networks

Recent years have witnessed the outstanding success of deep learning in various fields such as vision and natural language processing. This success is largely indebted to the massive size of deep learning models that is expected to increase unceasingly. This growth of the deep learning models is accompanied by issues related to their considerable energy consumption, both during the training and inference phases, as well as their scalability. Although a num-

ber of work based on unconventional physical systems, such as wave-based frameworks, addresses the issue of energy efficiency in the inference phase, efficient training of deep learning models has remained elusive. So far, training of digital deep learning models mainly relies on backpropagation, which is not suitable for physical implementation as it requires perfect knowledge of the computation performed in the so called forward pass of the neural network.

Here, we tackle this issue by proposing a simple deep neural network architecture augmented by a biologically plausible learning algorithm, referred to as 1

Introduction

Deep learning has emerged as a breakthrough technology with outstanding success in various fields such as vision, natural language processing (NLP), and speech recognition (1,2).

Although these algorithms attempt to mimic the functioning of the human brain, they are essentially executed on a software level using traditional von Neumann computing hardware. Nevertheless, artificial neural networks (ANNs) based on digital computing are currently experiencing challenges concerning energy consumption and processing speed (3). An example of the considerable energy usage involved in training language models can be seen in the case of GPT-3, which has 175 billion parameters and required 1.3 GWh of electricity during training (4). As a result, researchers are increasingly exploring the implementation of ANNs on alternative physical platforms, including but not limited to optics (5)(6)(7)(8)(9), spintronic (10,11), nanoelectronic devices (12)(13)(14), and photonic hardware (5,15), as well as acoustic systems (16,[START_REF] Momeni | Sixteenth International Congress on Artificial Materials for Novel Wave Phenomena (Metamaterials)[END_REF]. Currently, two primary methods are mainly used. The first requires the designing of hardware that implements trained mathematical transformations through strict operation-by-operation mathematical isomorphism, which is mainly aimed at the inference phase of deep learning [START_REF] Burr | [END_REF](19)(20)(21).

The second category, known as deep Physical Neural Networks (PNNs), involves training the physical transformations of the hardware directly to execute the desired computations. PNNs offer the potential to create more scalable, energy-efficient, and faster neural network hardware by leveraging physical transformations and foregoing the conventional software-hardware division (22,23).

In a similar fashion to digital deep learning, increasing the depth of physical networks (deep PNNs) holds great potential for enhancing performance as it leads to an exponential expansion in network expression capabilities (24). So far, the training of PNNs has predominantly relied on a method called backpropagation (BP), which has been highly successful in the training of the digital ANNs (25). Yet, there are several reasons why BP is not a suitable choice for PNNs, one of which is the complexity and lack of scalability in the physical implementations of BP operations in the hardware (26)(27)(28)(29). For example, in (26), the authors use the adjoint variable method to experimentally implement BP in photonic neural networks which is limited to lossless and reciprocal PNNs. Additionally, their method requires three rounds of light propagation, either from left to right or right to left, throughout the PNN, which remains unachievable in the majority of PNN setups. Most proposals for PNNs implement BP calculations on an external computer using a digital twin of the physical system, commonly referred to as in-silico training. This usually comes at the cost of speed and an increase in the energy consumption during training. Additionally, the model might not accurately represent the real physical sys- tem which can lead to a potential simulation-reality gap and unfaithful inference time prediction (6-8, 10, 13, 14, 20, 30, 31). Recent work has made some progress towards addressing these issues. The physics-aware training method based on BP (PA-BP) (22) is the current state-ofthe-art training framework that circumvents some problems of the in-silico methods. However, the Achilles heel of PA-BP is that it still requires a differentiable digital model for the backward pass. This imposes some limitations on the use cases of the PA-BP, such as the slow training speed, and the high power consumption due to the digital backward pass and the need for extra memory accommodating the backward model of the PNN on a digital computer. Furthermore, as becomes evident later, in the event that the physical system undergoes strong perturbations, reusing PA-BP-trained PNNs, or even their fine-tuned versions could be quite challenging and in some cases futile, necessitating training of the models from scratch.

Another significant drawback of BP is its reliance on having complete knowledge of the computations graph carried out during the forward pass to accurately compute derivatives (23,(32)(33)(34)(35). When a black box is inserted in the forward pass, BP becomes impossible. Therefore, researchers are seeking alternative training methods for PNNs. For example, an approach that has been recently explored for training physical networks is the augmented Direct Feedback Alignment (DFA) method (23), which aims to avoid the need for a differentiable digital model.

However, this method is only compatible with certain physical networks where the nonlinear and linear layers can be separated. Furthermore, determining the nonlinearity form for PNNs through optimization procedures is still an ongoing challenge.

Here, we propose a simple and physics-compatible PNNs architecture augmented by a biologically plausible learning algorithm, called the model-free forward-forward (MF-FF) training.

The proposed method enables direct training of arbitrary PNNs locally without the need to know the nonlinear physical layers and train a digital twin model. To the best of our knowledge this is the first work demonstrating a BP-free training of PNNs that circumvents all challenges as-sociated with the hardware-realization of BP. In this BP-free contrastive learning method, the standard backward pass, typically performed by a digital computer, is replaced with a single forward pass through a physical system. This substitution can significantly improve training speed compared to other hardware aware training framework, reduce digital computations and memory usage, and lower power consumption in the training phase of wave-based PNNs. In this work, we benchmark our learning framework against the state-of-the-art physics-aware training schemes proposed in the literature, i.e., in-silico and backpropagation methods. We show the robustness and adaptability of the proposed method compared to its sister schemes, even in systems exposed to unpredictable external perturbations. To showcase the universality of our approach, we perform experimental vowel and image classification using three wave-based systems that differ in terms of the underlying wave phenomenon and the type of non-linearity involved. Our first example consists of a chaotic acoustic cavity implemented with non-linear scatterers. The second example is a chaotic microwave cavity whose transfer function is massively parameterized by a programmable metasurface with structural non-linearity. Our third example is an optical multimodal fiber with readout non-linearity . Our approach results in high-accuracy achieving hierarchical classifiers that make use of the distinct physical transformations of each system and offers a pathway to significantly enhance the energy efficiency and speed of deep learning not only during inference but also during the training phase. each layer, is passed to the next layer. The subsequent layer then carries out the same process hierarchically on the output of the previous layer. This architecture shares some similarities with conventional deep Reservoir Computing (deep-RC) systems (36). However, here all augmented linear multiplications are trained, in contrast to the traditional deep-RC where only the final layer is trained (37,38). If nonlinear physical transformations possess some form of memory, the network can be categorized as a trainable deep-RC. Among the various contrastive learning approaches (35), the forward-forward algorithm (32), inspired by Boltzmann machines [START_REF] Hinton | Parallel distributed processing: Explorations in the microstructure of cognition[END_REF], was recently proposed by Geoffrey Hinton and has undergone further improvement through several proposals in computer science [START_REF] Zhao | [END_REF][START_REF] Ororbia | [END_REF]. Here, we implement the model-free version of the forward-forward algorithm, for the aforementioned PNNs-based architecture. Each nonlinear physical system performs a nonlinear transformation on input data, which can be expressed as

Model-free forward-forward training

h (l ) = f (l ) N (W (l ) p x (l ) ), where x (l ) , W (l )
p , and f (l ) N correspond to the physical inputs (e.g., optical intensity, electric voltage, vibration), physical interconnections (e.g., optical, electrical, or mechanical coupling) in the physical system, and physical nonlinearity (e.g., nonlinear optical, magnetic, or mechanical effects) in layer l , respectively. Indeed, W (l ) p and f (l ) N signify mixing operation and non-linear kernel of the l -th physical systems. Afterward, the output of layer l can be expressed as the multiplication of h (l ) by the augmented l ) . Such trainable matrix multiplications can be performed through either digital or physical systems, for instance using Mach-Zehnder Interferometer (MZI) integrated photonics (42) Spatial Light Modulators (SLMs) in optics (21,43). The goal here is to train W (l) t locally. Instead of a forward and backward pass, here, we use two physical forward passes: a positive and a negative forward pass through the physical systems, each running on different physical inputs. The positive physical pass, y

trainable weight matrix W (l ) t , y (l ) = W (l ) t h (
(l ) pos = W (l ) t f (l ) N (W (l ) p x (l ) pos ),
uses positive inputs that include the input dataset and correct labels, while the negative physical pass, y

(l ) neg = W (l ) t f (l ) N (W (l ) p x (l )
neg ), uses negative inputs that include the input dataset and incorrect labels (see Fig. 1a). In each layer, we calculate the so-called "goodness" function, defined as the sum of the squared activities for the positive and negative physical passes. For instance, in physical systems such as optics, the squared activities correspond to the optical intensity on the CCD camera. Eventually, for each layer l, W (l) t is trained by minimizing the following loss function

L (l ) = log 1 + exp -θ j y (l ) pos,j 2 - j y (l ) neg,j 2 (1)
where θ is a scale factor. During the inference phase, we input a particular label into the PNNs and accumulate the goodness values for all layers. This process is repeated for each label separately. The label with the highest accumulated goodness value is then selected as the output (see supplementary section S2 for more details). The proposed method is also capable of integrating non-differentiable physical systems or components between the layers.

Diverse PNNs for vowel and image classification

In Figure 1, we present three deep PNN classifiers for the vowel and Mnist (Modified National Institute of Standards and Technology database) handwritten digit classification task, based on three distinct physical systems. Although there have been proposals that explore wave-based analog computing for linear operations such as multiplication and convolution [START_REF] Zhou | [END_REF](44)(45)(46)(47)(48)(49)(50)(51), it is important to note that PNNs require nonlinearity to effectively handle a wide range of learning tasks. These tasks include regression and classification, which require nonlinear operations for accurate performance. We consider three wave-based physical systems, each of which features a unique origin of physical nonlinearity and underlying wave phenomenon, highlighting the diversity of physical networks that can be employed. We take advantage of the multiple scattering induced by the disordered environment in these physical systems (chaotic cavities and multimodal fibers) to leverage the mixing process. In addition to provide the required linear mixing, it also helps to intensify the overall nonlinearity (16).

Acoustic Chaotic Cavity with Non-Linear Scatterers

In acoustics, an air-filled multimode cavity composed of multiple nonlinear meta-scatterers randomly placed on the cavity top wall and multiple rigid scatterers inside the cavity is employed (materials and methods, supplementary section S3). The nonlinear meta-scatterers are designed based on an active control strategy. Herein, the positive and negative data are encoded onto the amplitude of each frequency component composing the excitation waveforms, that are then injected into the nonlinear system through loudspeakers positioned on the right side of the cavity. The output of the physical system is measured using microphones below the metascatterers.

We investigate the vowel classification performance of two layers acoustic-PNN (see Fig. 2a).

In order to compare the results of MF-FF with ideal BP, and in-silico training, we accurately model the forward pass of acoustic-PPN by a digital neural network (see the supplementary section S3). When the acoustic-PNN is trained using MF-FF, it performs the classification task with 97.62% and 94.23% train and test accuracy, respectively (see Figs. 2d-f). Figure 2c shows the comparison of the classification results obtained for MF-FF, ideal BP, and in-silico training.

A schematic visual representation of the aforementioned methods is provided above Fig. 2c.

The complete comparison between different hardware-based training methods is provided in supplementary section S2.

As evidenced by Fig. 2c, in-silico training performs poorly, reaching only a maximum vowel classification accuracy of ∼60%. When there is a gap between the reality and the simulation of a physical system ( called the reality-simulation gap), the accuracy of inference will decrease. In contrast, MF-FF succeeds in accurately training the acoustic-PNN, performing similarly to the ideal BP algorithm used as a baseline. The key advantage of MF-FF stems from the execution of both forward passes through the physical hardware, rather than simulations.

Microwave Massively Parametrized Chaotic Cavity with Structural Non-Linearity

In the microwave regime, we leverage a "structural non-linearity" such that we can implement nonlinear mathematical operations at low power levels with a linear scattering system. Our starting point is an irregularly shaped electrically large metallic enclosure with strong modal overlap. It is coupled via two coax-to-waveguide adapters to two asymptotic scattering channels, and the transfer function between these two channels can be measured using standard equipment such as a vector network analyzer. We then massively parametrize this cavity by covering one of its walls with a programmable metasurface. For each meta-atom and each polarization, we can individually configure the local boundary conditions of the cavity with 1-bit precision (two possible states). Our setup is shown in Fig. 3 was recently used to implement with high fidelity and in situ reprogrammability desired linear transfer functions for signal differentiation (52) and routing (53). In the present work, however, we seek a non-linear mapping. Hence, we define the metasurface configuration as the input and the transfer function as the output of our mathematical operation. Indeed, this relation is in general nonlinear due to the mutual coupling between meta-atoms caused by their proximity and, more importantly, the reverberation (54). While previous work in Ref. (49) sought to limit the reverberation in order to implement a linear transformation with the metasurface configuration as input and the transfer function as output, here we deliberately seek to maximize the reverberation to boost the nonlinearity (further explained in the supplementary section S4). Incidentally, this type of reverberation-induced structural nonlinearity was recently also transposed to the optical domain (55). 3c-e).

Optical Multimode Fiber with Readout Non-Linearity

In the optics part, we use an optical system that comprises an SLM, a scattering medium consisting of a step-index multimode fiber (MMF), and a CMOS camera (see Fig. 1d, materials and methods). In this setup, the positive and negative datasets are encoded onto the SLM, and after passing through the MMF, the resulting optical intensity is measured on the camera. The physical optical system performs a complex spatiotemporal transformation. This transformation involves the propagation of spatially modulated laser pulses through an MMF. The propagation of an ultrashort pulse inside a MMF is a highly complex process that involves spatial and temporal interactions of electromagnetic waves coupled to hundreds of different propagation modes (56,57). Although this transformation is linear in the complex domain, the process becomes nonlinear due to the data being encoded onto the phase (SLM) and the subsequent measurement of the intensity squared on the camera. Here, we use a two-layer optics-PNN to perform classification tasks on two different datasets: Vowel and Mnist (see Fig. 22) and the proposed MF-FF method is presented, with the focus on their ability to recover the classification accuracy after applying perturbation.

digital processor to a fast, energy-efficient hybrid physical-digital processor, showcasing the potential for optimizing both performance and energy efficiency in machine learning applications.

After posting the initial version of this work on arXiv, we became aware that a related idea and experiment showed for optics was being pursued by another team, see preprint at (58).

Real-Time Adaptable Learning

We now aim to show the superior robustness of MF-FF in the context of real-time and adaptable learning, where the physical data transformer may undergo changes due to slow dynamics of the physical system during the runtime or external hard perturbations. Let us consider a deep optics-PNN with six layers, as depicted in Fig. S6a, which has already been trained on vowel datasets and is currently in the inference phase. The transformation function of each physical system is f 0 (θ), where θ is the physical input. We perturb the physical systems at a specific time (examples of such perturbations include changes in the MMF state or the positions of lenses or masks, etc), which results in a change in the transformation function of each physical system from f 0 (θ) to f p (θ) (see Fig. S6b). We simply perturb the transmission matrix of the optical setup by adding Gaussian noise with mean µ and standard deviation σ to simulate this situation.

As we observe in Fig. S6d, the test accuracy drops as expected after applying the perturbation.

The question now is whether the training method can restore the accuracy by retraining the optics-PNN after some epochs. We compare our results with the PA-BP method (22) which uses a digital model for the backward pass and the physical system for the forward pass. As shown in Fig. S6c, PA-BP struggles to restore accuracy with increasing perturbation intensity. For instance, the test accuracy oscillates around 55% for a small perturbation (red dots in Fig. S6c)

and worsens further for more intense perturbations. In stark contrast, the proposed MF-FF can easily recover accuracy after a few epochs, regardless of the intensity of the perturbation applied. (see Fig. S6d).

This striking adaptability owed to the fact that MF-FF executes both forward passes through the physical hardware, rather than digital models. In contrast, the PA-BP method uses a digital model that becomes completely inaccurate after hard-perturbation, necessitating re-training with a new dataset and wasting energy. This study highlights the robustness and adaptability of the proposed model, even in dynamic or unpredictable environments. Additionally, one prominent advantage of MF-FF training is its robustness to an external perturbation. We observed that the MF-FF training is superior not only to the in-silico training but also to the current state-of-the-art PA-BP method. The performance gap between our method and the PA-BP continues to widen as the depth of the NN increases. This is because in PA-BP training, the mismatch between the gradients of the non-perturbed and perturbed systems accumulate over the depth of the ANN. Therefore, the more layers a network has, the more sensitive it is to an external perturbation. This is in contrast to the MF-FF trained framework in which each layer of the network is trained separately without using a model or direct dependence to other layers that is required for BP.

Discussion

Reference (21) has recently revealed the exceptional efficiency of optical transformer models during the "inference" phase (achieving a > 8, 000× energy-efficiency advantage over state-of-the-art digital-electronic processors ). The training method proposed in this paper represents a viable candidate for "training" such optical transformer models, enabling the potential to achieve energy-efficiency and speed advantages. We further study the potential of PNNs in terms of reducing energy consumption and improving the computing speed compared to digital training/inference schemes of ANNs in Section S6 of the supplementary material.

We conclude that even though PNNs have the potential to improve deep neural networks' speed and energy consumption to a large extent, a number of challenges still need to be addressed before PNNs can replace their digital counterparts. For example, it is not yet clear how some widely used mathematical operations, such as normalization units, including layer/batch normalization, could be implemented in hardware. Second, how the physical systems could be compactly scaled to billions of parameters. Therefore, we can expect, at least in short term, that PNNs will enable efficient hybrid combinations of in-hardware training/inference rather than a complete replacement for digital processing units.

scattering channels. The 2-bit programmable metasurface contains 76 meta-atoms that cover 8 % of the cavity surface and efficiently modulate the field inside the cavity within a 400 MHz interval centered on 5.2 GHz. The two bits of control per meta-atom are assigned to the two orthogonal field polarizations. Using a vector network analyzer (Rhode & Schwarz ZVA 67), the transmission spectrum between the two ports is measured. The cavity has a composite quality factor on the order of 370 and around 23 modes overlap at a given frequency (52,53). The system's scattering matrix is subunitary due to significant Ohmic losses on the metallic walls. We randomly assign the available 152 1-bit programmable degrees of freedom of our programmable metasurface to 40 groups, coined macro-pixels, since our input vector's dimensions are 1 × 40.

All meta-atoms within a given macro-pixel are configured identically. We select twenty decorrelated frequencies within the metasurface's operation bandwidth and consider the transfer function intensity at these frequencies as outputs. The relation between inputs and outputs is hence non-linear due to the structural non-linearity (54) and the readout non-linearity. Our use of this microwave setup is distinct from previous works on over-the-air microwave analog computing with programmable metasurfaces in chaotic cavities. While, on the one hand, Ref. (49) also defined a mathematical operation in which the inputs were related to the metasurface configuration and the outputs to the transfer function, it aimed at implementing linear operations and hence sought to minimize the reverberation-induced structural non-linearity whereas we seek to maximize it in the present work. Refs. (52,53), on the other hand, considered the system's transfer function as the mathematical operation in order to implement desired linear transfer functions for high-fidelity in situ programmable signal differentiation and routing; the inputs in Refs. (52,53) were hence the incident wavefronts rather than the metasurface configuration.

gradient-free optimization techniques can be used to update parameters in each layer, it may impose some limitations, such as slow convergence rates, getting stuck in local minima, and requiring heavy digital computations.

Algorithm 1

Input: Training dataset (x i , y i ), N exter and N inter (number of external and internal training epochs, respectively), L (number of layers), θ (scale factor), and PNN forward (forward pass of PNN) .

Output:

← W l h pos , W l h neg L ← log 1 + exp -θ j h (l) pos,j 2 -j h (l) neg,j 2
Updating W l by minimizing L end for

To mitigate such limitations, we employed a simple deep PNN architecture as illustrated in Fig. 1a the measured front pressure and according to the following control law

i(t) = Bl S d G NL |p f (t)| α NL , (2)
where G NL and α NL are two tunable parameters that produce nonlinearity in the system, and

Bl and S d are the force factor and effective cross-sectional area of the speaker, respectively.

The optimal values for high nonlinearity and control stability for G NL and α NL are reported in Table S1.

The measurement is performed 10,000 times with different amplitude weightings for each run and each input. The set of input and output data is thus composed of two matrices with 40 × 10000 elements.

S3-2. Input-output transformation characterization

To compare our results from the MF-FF method with other training methods such as ideal BP and in-silico training methods, we trained a digital neural network for the forward pass of the acoustic system (see Fig. 2c of the main paper). We collected training data by recording the outputs of the physical system for a selection of N input vectors, here N was 10 4 and the input vectors were drawn from a uniform distribution. The data set was divided into training and validation sets and used to train the neural network. A fully-connected neural network was used to map the input-vector with the dimension of 1 × 20 to the output-vector with the dimension of 1×20. The numbers of hidden layers are three with 100, 200, and 100 neurons, respectively. We used the Sigmoid Linear Unit (SiLU) activation function and Adam optimizer. After training allowing them to maintain high classification performance even in the presence of stochastic failures in any individual sub-component (61).

In order to implement the results of in-silico (see Fig. 2c of the main paper), we added small noise (Gaussian noise with a mean of zero and standard deviation of 0.025) to the parameters of the forward model of the nonlinear acoustic system. As long as there is a small gap between the reality and the digital model of the physical system, the accuracy of inference will decrease, as illustrated in Fig. 2c of the main paper. puts is hence non-linear due to the structural non-linearity (54) and the readout non-linearity.

Our use of this microwave setup is distinct from previous works on over-the-air microwave analog computing with programmable metasurfaces in chaotic cavities. While, on the one hand, Ref. ( 49) also defined a mathematical operation in which the inputs were related to the metasurface configuration and the outputs to the transfer function, it aimed at implementing linear operations and hence sought to minimize the reverberation-induced structural non-linearity whereas we seek to maximize it in the present work. Refs. (52,53), on the other hand, considered the system's transfer function as the mathematical operation in order to implement desired linear transfer functions for high-fidelity in situ programmable signal differentiation and routing; the inputs in Refs. (52,53) were hence the incident wavefronts rather than the metasurface These values evidence that there are significant non-linear effects in the mapping from metasurface configuration to transmission spectrum (which we refer to as structural non-linearity) that were shown in Ref. (54) to originate from two factors: (i) mutual coupling due to proximity between close-by meta-atoms; and (ii) reverberation-induced long-range coupling between all meta-atoms. Again, we have highlighted in Fig. S5b the 20 frequencies that we selected for our analysis in the main text.

S4-2. Input-output transformation characterization

Similar to the procedure outlined for the acoustics-PNN, we trained a digital neural network to serve as a surrogate forward model that approximates the forward pass of our microwave system (see Figure 3a of the main paper). We obtained training data by recording the physical Section S6. Energy consumption and computing rate analysis At the inference stage, the energy consumption can be broken down into the energy cost of the data loading into the physical system, a feedforward pass of the wave fields in the system, and a readout phase into the digital system (21,26). Therefore,

E inference = E load + E feedforward + E readout (3)
The E load consists of the energy required to read from the memory, digital to analog (DAC) conversion as well as the modulation energy of the wave modulator (E load = E read + E DAC + E mod ). Similarly, the E readout is comprised of a wave amplification, analog to digital con-version (ADC), and writing to memory (or energy consumption of camera). The E feedforward consists of the energy of the waves in the system including the wave sources in addition to maintaining any active components in the system including SLM or DMD. The details of energy consumption are listed in Table S2 for the proposed optical system.

At the training phase, the energy consumption consists of the energy at the inference E inference plus the energy consumption of a single gradient update of the learnable parameters. The E gradient consists of a digital gradient update to the learnable parameters of the physical network as well as an additional E load to implement a physical update to the parameters. Therefore,

E training = 2E load + E feedforward + E readout + E gradient (4)
In the above analysis, the energy contributions of the E load , E readout and E feedforward scale linearly with the number of the data points N (image's pixels, vowel's frequencies, etc.) loaded/read out of the system (we assume the system has the same number of inputs and outputs) while E gradient having an N 2 scaling contribution.

The energy consumption for the training/inference of a fully digital neural network is measured in the total number of digital floating-point operations per second (Flops) or the total number of multiply-accumulate (MAC) operations of a given vector-matrix product. In our proposed FF-model, the digital energy cost of a forward pass through the linear layer requires quadratic N 2 MAC operations while a PNN's energy scales linearly, i.e. N(E load + E readout ).

While the cost of E feedforward is virtually free, E load , E readout are the leading contributors to the total energy consumption of a PNN. However, in large models, these extra costs amortize the per-MAC cost of the PNN compared to the constant per-MAC cost of a digital model.

While PNNs hold great potential for decreasing the energy cost of the computation, the computing rate of a PNN could be unpleasantly slow. Wave modulators are usually bottlenecking this information processing speed. Spatial light modulators (SLMs) and digital micromirror devices (DMDs) are currently the major optical modulators with speeds on the orders of Hzs and kHzs. For inference, the effective computation rate is limited by the rate of the input data fed to the system (∼ GHz) whereas for training, the rate is limited to the refresh rates of modulators containing the learnable weights of the system. Switching

Figure 1 -

 1 Figure 1 -Deep physical neural networks. a, A simple and physics-compatible deep neural network that employs a sequence of nonlinear physical data transformers augmented by trainable matrix multiplications, trained by the model-free forward-forward (MF-FF) technique. At each layer, the nonlinear physical data transformer executes nonlinear mapping between the input and output space, and the MF-FF trains the trainable linear multiplications to determine the optimal decision boundary for positive and negative data. We consider three physical systems that vary in terms of the underlying wave phenomenon and the type of non-linearity. b, In acoustics, input data is encoded into the intensity of sound waves at different frequencies injected on the left side of the cavity. Sound waves propagate through a chaotic cavity that comprises multiple rigid cylindrical diffusers and nonlinear membranes. The transformed waveforms are received by multiple microphones. c, In the chaotic microwave cavity, input data is encoded into the programmable metasurface configuration inside the metallic disordered cavity. The outputs are obtained from the waves' spectra (transfer function). d, In optical setup, input data is encoded onto the SLM, and after passing through the MMF, the resulting optical intensity is measured on the CCD camera.

Figure 1a shows aFigure 2 -

 2 Figure1ashows a simple and physics-compatible deep PNN including N nonlinear physical data transformers augmented by trainable linear multiplications. Each nonlinear physical data transformer performs a nonlinear mapping between the input and output space "effortlessly", followed by the use of augmented trainable linear multiplication to classify distinct classes via a local and contrastive training algorithm. The orientation of neuron activities, the output of

Figure 3 -

 3 Figure 3 -Microwave-PNN. a, The topology of microwave-PNN consists of a three-layer PNN with skip connections. Each layer comprises a microwave-PNN augmented by trainable matrix multiplication. b, Photograph of the experimental setup. c, The train and test classification accuracy versus training epoch for the vowel recognition task. d and e, The confusion matrix for the PNN on the train and test set, respectively.

Figure 4 -

 4 Figure 4 -Optics-PNN. a, The topology of optics-PNN consists of a two-layer PNN. Each layer comprises an optics-PNN augmented by trainable matrix multiplication. b and f, An example of input data, including raw data and its label representation on SLM, along with its corresponding output on a CCD camera for vowel and Mnist datasets, respectively. c and g, The train and test classification accuracy versus training epoch for the vowel recognition and Mnist tasks, respectively. d and e, The confusion matrix for the PNN on the train and test set for the vowel recognition task, respectively. h and i, similar results for the Mnist task.

Figure 5 -

 5 Figure 5 -Robustness of deep PNNs against unpredictable external perturbations. a, A deep PNN consists of six layers of optics-PNN augmented by trainable matrix multiplication. The deep PNNs trained on vowel datasets and is currently in the inference phase b, Applying hard-perturbation by adding Gaussian noise with the mean of µ and standard deviation of σ to the transmission matrix of MMF. c and d, A comparison between PA-BP (22) and the proposed MF-FF method is presented, with the focus on their ability to recover the classification accuracy after applying perturbation.

  Due to the unprecedented growth in the size of ANNs, including models such asGPT 4 with around hundreds of billion parameters, the cost of both the training and inference phases of these networks has increased exponentially. Training of such massive ANNs is particularly extremely expensive and time-consuming, with training times taking up to several months. Specialized hardware such as PNNs have the potential to drastically decrease these costs. A few methods have been proposed for training PNNs, but they all face issues such as mismatch between the forward model and the physical system or robustness issues. This is because these methods perform the entire backward pass through a digital computer during training, involving either a digital model in PA-BP or numerical simulations in in-silico training, which can hinder their effectiveness in the training phase. MF-FF enables forward passes through physical systems, resulting in a significant speed-up not only during inference but also in the training phase.

Figure S2 -

 S2 Figure S2 -Active control scheme and set-up photography a, Photography of the experimental set-up with the localization of the active nonlinear resonators. b, nonlinear active control scheme. c, Closed-up photography on the active nonlinear membrane (ANM).

Figure S4 -

 S4 Figure S4 -Details about the programmable metasurface. a, Close-up view. b, Setup to characterize normal-incidence scattering properties of the metasurface. c, Normal-incidence characterization of the scattering of the two different meta-atom configurations.

FrequencyFigure S5 -

 S5 Figure S5 -Additional characterization of the microwave setup. a, Field stirring efficiency. b, Linearity metric.

Figure S6 -

 S6 Figure S6 -Comparison of the forward model and physical microwave-PNN for four randomly selected inputs from the test dataset.

Figure S7 -

 S7 Figure S7 -Transformation of optics-PNN for ten randomly selected digits from the Mnist dataset.

  from an LC SLM to a DMD could improve the inference speed by nearly 3 orders of magnitude, reaching up to 23000 frames/second. Although the slow computation rate of the training is amortized by the fast rate of the inference (training is done only once), improvements to electronics of the modulators could further improve the training efficiency of PNNs in the future.

  Trainable parameters of matrix W l for l th layer Initialization Parameters W l for Epoch ∈ {1, 2, 3, ...N exter } do Creating x pos and x neg from (x i , y i ) h pos , x pos ← h neg , x neg for layer ∈ {1, 2, 3, ...L} do h pos , h neg ← PNN forward (h pos ), PNN forward (h neg ) h pos , h neg ← TrainLayer(h pos , h neg ) end for end for Func: TrainLayer (h pos , h neg ) for Epoch ∈ {1, 2, 3, ...N inter } do h pos , h neg

Table S1 -

 S1 Control & speaker parameters.

	G NL α NL Bl (Tm) S d (m 3 )
	ANM1 5.7	1.6	1.1	0.0012
	ANM2 5.5	1.7	1.0	0.0012
	ANM3 5.3	1.5	1.0	0.0012
	ANM4 5.6	1.4	1.1	0.0012

Table S2 -

 S2 The energy budget breakdown of the optical computing system. stages separately. The FF training is also particularly suitable for breaking down the energy estimation per layer.

		Energy Consumption	Details and/or Ref
	E read	1 pJ/bit for off-chip memory,	(21, 63)
		and 0.3 pJ/bit for SRAM
	E DAC	10 pJ per 5-bit sample at 10	This is achievable with 100
		GHz	mW at 30.1dB SFDR (21)
	E mod	45× r t pJ /MAC	by considering the 1mW
			continuous laser beam. r t is
			the refresh time correspond-
			ing to the refresh time of
			SLM or DMD. (59)
	E SLM	15nJ /MAC	Assuming 2W for operation
			of a 1920×1152 SLM, work-
			ing at 60 Hz (21)
	E DMD	94pJ /MAC	For	Texas	Instruments
			DLP9500 DMD model (with
			a 1920 × 1080 resolution and
			23148 patterns per second
			) (57)
	E ADC	3.17 pJ per 7-bit sample	(21)
	E readout	1.4 nJ /MAC	For LUX 1310 camera (57)

Here we discuss the energy cost of the physical neural network proposed in this manuscript. We analyze the energy cost of a layer of the physical neural network for the inference and training

Acknowledgments

A.M. and R.F. acknowledge funding from the Swiss National Science Foundation under the Eccellenza grant number 181232. P.d.H. and R.F. acknowledge funding from the ANR-SNF PRCI program (project "MetaLearn"; ANR-22-CE93-0010-01). B.R. acknowledges that all materials pertaining to the optical experiment presented in this manuscript have been sourced from previously published material that is publicly available and was acquired during his tenure at the Laboratory of Applied Photonics Devices at EPFL. This work has no connection to his current employer in any capacity or form. A.M. and M.M. thank Dr. X. Guo for useful discussions.

Contributions

A.M. conceived the idea, designed the computational engine, and carried out both the theoretical and numerical simulations as well as a part of the acoustic experiment. B.R. provided the optics data and interpretation of machine learning results. M.M. carried out the acoustic experiment. P.d.H. carried out the microwave experiment. R.F. supervised the project. All authors contributed to the interpretation of the results and the writing of the manuscript.

Supplementary materials

Materials and Methods

Supplementary Text

Figs. S1 to S7 Each delivers a different linear waveform with a 10 frequencies content in the 350 -500 Hz range with randomly weighted amplitudes. The output of the system consists of 4 highly nonlinear pressure waveforms measured at the output of the active control at 4 different locations.

The measurement is performed 10,000 times with different amplitude weighting for each run and each input.

S1-2. Microwave system

Our microwave system is shown in Fig. 3b of the main paper and consists of an irregularly shaped electrically large metallic cavity (0.385 m × 0.422 m × 0.405 m; 0.0658 m 3 ) whose scattering properties are massively parametrized through a programmable metasurface. Two waveguide-to-coax adapters (RA13PBZ012-BSMA-F) connect the cavity to two asymptotic

S1-3. Optical system

The data used for the optical experiment are from the published dataset in Rahmani et al. (59).

The optical system therein consists of a spatial light modulator, a scattering medium of a stepindex multimode fiber of length 0.75m and core diameter 50µm with an aperture size 0.22, and a CMOS camera. The entire system is operating at the 532nm and 1 mW power continuouswavelength light source corresponding to ∼ 1050 modes of the fiber for one polarization. Using interferometry measurements (60), a transmission matrix of the optical system is obtained which allows mapping an input 2-dimensional optical field to its complex output field. The columns of the transmission matrix contain the response functions of the system for each of the modes of the fiber, allowing to faithfully calculate the optical responses of the system to arbitrary inputs.

Section S2. Comparison of different PNN training algorithms

In Figure S1, a comprehensive comparison of different algorithms including ideal BP, in-silico, PA-BP, and MF-FF is depicted.

All algorithms, except MF-FF, utilize backpropagation as the core technique and consist of four key steps: forward pass, error vector computation, backward pass, and parameter update.

The key distinction among these algorithms lies in the choice of physical transformation, either

for the forward and backward passes. For example, the ideal BP uses the exact function of the system f exact for both the forward and backward passes. However, in-silico employs f model for both the forward and backward passes. Physics-aware BP utilizes f physical for the forward pass and f model for the backward pass.

In contrast to the aforementioned methods for training PNN, MF-FF employs two forward passes without employing any backpropagation through layers. The original paper was recently proposed by Geoffrey Hinton (32). The concept involves substituting the traditional forward and backward passes of backpropagation with two forward passes that function in the same manner, but with contrasting objectives and operate on distinct datasets. The positive pass involves processing real data and adjusts the weights to increase the goodness in every hidden layer. On the other hand, the negative pass involves processing negative data and adjusting the weights to decrease the goodness of each hidden layer. In this paper, we use the balanced contrastive loss which is recently proposed in Ref. (34) (Eq. 1 in the main paper). However, other forms of imbalance positive-negative losses are possible (32).

To update the parameters of each layer in conventional fully-connected multilayer neural for 200 epochs, the model achieved a mean squared error of 2.8 × 10 -3 on the test dataset. In Fig. S3, we randomly selected four samples from the testing dataset and compared the model's predictions to the ground truth physical output. This visual comparison allows us to assess how well the model predict the output of physical system.

The proposed acoustics-PNN consists of two layers with skip-connections (see Fig. 2a of the main paper) and was trained by MF-FF. We found that incorporating skip connections into the network was beneficial. This feature is particularly beneficial as the acoustic transformation lacks support for identity operation, making it susceptible to information loss. These connections, inspired by residual neural networks, were added to allow the network to act further as an ensemble of sub-networks (22,61), thus this contributes to the robustness of the architectures, Section S4. Microwave-PNN S4-1. Experimental setup

Our microwave system is shown in Fig. 3b of the main text and consists of an irregularly shaped electrically large metallic cavity (0.385 m×0.422 m×0.405 m; 0.0658 m 3 ). Two waveguide-tocoax adapters (RA13PBZ012-BSMA-F) connect the cavity to two asymptotic scattering channels. Using a vector network analyzer (Rhode & Schwarz ZVA 67), the transmission spectrum between the two ports is measured. The cavity has a composite quality factor on the order of 370 and around 23 modes overlap at a given frequency (52,53). The system's scattering matrix is subunitary due to significant Ohmic losses on the metallic walls.

The cavity's scattering properties are massively parametrized through a programmable metasurface. A close-up view of the programmable metasurface (purchased from Greenerwave) is shown in Fig. S4a. The 2-bit programmable metasurface contains 76 meta-atoms that cover 8 % of the cavity surface and efficiently modulate the field inside the cavity within a 400 MHz interval centered on 5.2 GHz. The two bits of control per meta-atom are assigned to the two orthogonal field polarizations ("H-Pol" and "V-Pol" in the inset in Fig. S4a). The working principle of the 1-bit control over each polarization follows that proposed in Ref. (62). At the central working frequency (around 5.2 GHz), the metasurface can tune the reflected phase by roughly π, as seen in Fig. S4c, implying that it can roughly mimic Neumann or Dirichlet boundary conditions.

We randomly assign the available 152 1-bit programmable degrees of freedom of our programmable metasurface to 40 groups, coined macro-pixels, since our input vector's dimensions are 1 × 40. All meta-atoms within a given macro-pixel are configured identically. We select twenty decorrelated frequencies within the metasurface's operation bandwidth and consider the transfer function intensity at these frequencies as outputs. The relation between inputs and out-configuration.

Finally, we now detail two additional characterizations of our setup.

First, in Fig. S5a, we quantify the control of the programmable metasurface over the transmission spectrum between our two antennas in terms of the K-factor that is commonly used in electromagnetic compatibility and wireless communications. The K-factor is defined at each frequency as the power ratio of the unstirred to the stirred field components:

, where µ(f) and σ(f) are the mean and standard deviations of an ensemble of transmission spectra corresponding to random metasurface configurations. The lower the Kfactor is, the more efficiently the metasurface stirrs the field, i.e., the more control it has over the transmission between the two antennas. We have highlighted in Fig. S5a the 20 frequency (regularly spaced) that we selected for our analysis in the main text. The stirring efficiency is seen to be quite frequency dependent, since it is determined both by the metasurface's frequency response and by the frequency-dependent scattering properties of the cavity.

Second, in Fig. S5b, we quantify the linearity of the complex-valued transmission spectrum's dependence on the metasurface configuration. Recall that in the main text we only work with the intensity rather than the complex-valued fields analyzed here, which adds an additional readout non-linearity. Here, we seek to only quantify the structural non-linearity. We use the linearity metric ζ introduced in Ref. (54). To compute this metric at a given frequency, we begin by fitting the best possible linear model, h = h 0 + t T c, to the data set of random metasurface configurations c and the corresponding transmission coefficients h. Then, we compute the linear model's error on unseen data, ∆h, and evaluate an SNR-like metric that treats this prediction error like noise in the usual SNR definition: ζ = SD i (h i )/SD i (∆h i The higher the value of the linearity metric ζ is, the more linear is the relation between the transmission spectrum and the metasurface configuration. We observe that the value of ζ is lowest in the middle of the metasurface's frequency interval of operation, and that it is always in the vicinity of 10 dB.

The proposed microwave-PNN consists of three layers with skip-connections (see Fig. 3a of the main paper) and was trained by MF-FF.

Section S5. Optics-PNN S5.1 Experimental setup

The experimental system from which the optical dataset is obtained is explained in detail in the previous work (59).

S5-2. Input-output transformation characterization

The MMF transformation consists of three steps: first, a Fourier transform is applied to the input image; second, the result is multiplied by a transmission matrix; and third, an inverse Fourier transform is applied to the product. Figure S7 displays the transformation of optics-PNN for ten randomly selected digits from the Mnist dataset.

The proposed optics-PNN consists of two layers (see Fig. 4a of the main paper) and was trained by MF-FF. The dimension of matrix multiplications is 51 × 51 and 676 × 676 for the Vowel and Mnist datasets, respectively.


domain: phys spi

text:
Malte Röntgen

Olivier Richoux

Georgios Theocharis

Christian V Morfonios

Peter Schmelcher

Del Philipp

Hougne

Vassos Achilleos

Equireflectionality and customized unbalanced coherent perfect absorption in asymmetric waveguide networks

We explore the scattering of waves in designed asymmetric one-dimensional waveguide networks. We show that the reflection between two ports of an asymmetric network can be identical over a broad frequency range, as if the network was mirror-symmetric, under the condition of so-called latent symmetry between the ports. This broadband equireflectionality is validated numerically for acoustic waveguides and experimentally through measurements on microwave transmission-line networks. In addition, introducing a generalization of latent symmetry, we study the properties of an N -port scattering matrix S. When the powers of S fulfill certain relations, which we coin scaled cospectrality, the setup is guaranteed to possess at least one zero eigenvalue of S, so that the setup features coherent perfect absorption. More importantly, scaled cospectrality introduces a scaling factor which controls the asymmetry of the incoming wave to be absorbed. Our findings introduce a novel approach for designing tunable wave manipulation devices in asymmetric setups. As evidenced by our acoustic simulations and microwave experiments, the generality of our approach extends its potential applications to a wide range of physical systems.

I. INTRODUCTION

Scattering of waves is unambiguously of fundamental importance in physics, finding applications in fields as diverse as high-energy physics [START_REF] Peskin | An Introduction to Quantum Field Theory[END_REF], X-ray diffraction [START_REF] Cullity | Elements of X-Ray Diffraction[END_REF], wave localization [START_REF]Introduction to Wave Scattering, Localization and Mesoscopic Phenomena[END_REF] or wave filtering [START_REF] Macleod | Thin-Film Optical Filters[END_REF][START_REF] Dokumacı | Duct Acoustics: Fundamentals and Applications to Mufflers and Silencers[END_REF][START_REF]Springer Handbook of Acoustics[END_REF]. In any scattering problem, the system's symmetries have a strong influence. For a reflection symmetry-the simplest geometric symmetry-, for instance, waves sent into the system from the two opposite sides of the symmetry axis (or plane) act identically, thus the two reflection coefficents are strictly equal.

Here, on the other hand, we are interested in exploring the scattering properties of geometrically asymmetric systems. In particular, we focus on systems featuring the recently introduced latent symmetry [START_REF] Smith | Hidden symmetries in real and theoretical networks[END_REF][START_REF] Kempton | Characterizing cospectral vertices via isospectral reduction[END_REF][START_REF] Röntgen | Latent Symmetries: An Introduction, Other[END_REF]. Such a symmetry is usually not visible in the original setup, but it becomes apparent after a suitable dimensional reduction, the so-called isospectral reduction [START_REF] Bunimovich | Isospectral Transformations: A New Approach to Analyzing Multidimensional Systems and Networks[END_REF]. What is interesting about a latent symmetry is its strong impact on the eigenmodes of the underlying system [START_REF] Röntgen | Latent symmetry induced degeneracies[END_REF]. A latent reflection symmetry, for instance, induces local parity on the eigenvectors of the underlying matrix M describing the system (for instance, the Hamiltonian or the scattering matrix) [START_REF] Kempton | Characterizing cospectral vertices via isospectral reduction[END_REF][START_REF] Morfonios | Cospectrality preserving graph modifications and eigenvector properties via walk equivalence of vertices[END_REF][START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF].

Interestingly, the impact of a latent symmetry goes beyond the system's eigenvectors and manifests itself also in certain relations of the powers of the matrix M [START_REF] Röntgen | On symmetries of a matrix and its isospectral reduction[END_REF]. For a latent reflection symmetry with respect to two sites u, v, the corresponding diagonal elements M k u,u and M k v,v of the matrix powers of M are the same for all positive integers k [START_REF] Kempton | Characterizing cospectral vertices via isospectral reduction[END_REF]. On a fundamental level, this result is interesting, as it shows a deep connection between the powers of M and its eigenvectors [START_REF] Morfonios | Cospectrality preserving graph modifications and eigenvector properties via walk equivalence of vertices[END_REF].

In the first part of this work, we carry the concepts of latent symmetry and of matrix powers to the realm of wave scattering. We start by applying the concept of latent symmetry to design geometrically asymmetric systems whose scattering properties carry the same traits as a reflection symmetric system, that is, broadband equireflectionality. We validate this intriguing property numerically in acoustic waveguide networks and experimentally in microwave transmission-line networks.

Then, in the second part of this work, we explore the impact of relations in the matrix powers of a general scattering matrix. Specifically, we show that an N -port system can be designed to feature coherent perfect absorption (CPA) using scaled matrix power relations. CPA implies the complete and irreversible transduction of the incident wave energy into other degrees of freedom such as heat, which is possible whenever the scattering matrix has a zero eigenvalue and the corresponding eigenvector is used as incident wavefront [START_REF] Chong | Coherent perfect absorbers: Time-reversed lasers[END_REF]. This condition is valid irrespective of the complexity of the wave system, applying to simple [START_REF] Landy | Perfect metamaterial absorber[END_REF] or disordered [START_REF] Imani | Perfect absorption in a disordered medium with programmable meta-atom inclusions[END_REF] systems excited by a single channel as well as simple [START_REF] Wan | Time-reversed lasing and interferometric control of absorption[END_REF] or disordered [START_REF] Pichler | Random antilasing through coherent perfect absorption in a disordered medium[END_REF] systems excited by multiple channels. In asymmetric disordered systems, the necessary wavefront is generally very complex unless the system is optimized to impose CPA with a specific wavefront [START_REF] Del Hougne | Coherent wave control in complex media with arbitrary wavefronts[END_REF]. Here, we show that a 3-port system designed to feature scaled cospectral-ity will enable CPA with a prescribed imbalance of the corresponding wavefront that is controlled by tuning the associated scaling factor.

This paper is organized is follows. Section II is dedicated to the concept of equireflectional scattering, including theoretical and experimental verifications of the theory. Section III focuses on the powers of the scattering matrix of a generic N -port setup. In particular, we introduce scaled matrix power relations which can be directly related to CPA. We then focus on a simple three-port system, which we optimize to feature scaled cospectrality and hence CPA with custom imbalance of the CPA wavefront. Finally, we conclude our work in Section IV.

II. SCATTERING OFF LATENTLY SYMMETRIC NETWORKS

In the following, we will discuss how a latent reflection symmetry leads to equireflectional scattering. After introducing this phenomenon by means of a simple example in Section II A, we show in Section II B how this phenomenon can be explained and designed; in Section II C and Section II D we numerically and experimentally validate our findings with acoustic waveguide and microwave transmission-line networks, respectively.

A. Setup and a first example Consider a general reciprocal two port system, with only a single mode traveling within each lead/port as depicted in Fig. 1(a). This problem is characterized by the following scattering matrix equation

ψ - 1 ψ - 2 = r 1 t t r 2 ψ + 1 ψ + 2 = S ψ + 1 ψ + 2 , (1)
where (ψ - 1 , ψ - 2 ) T and (ψ + 1 , ψ + 2 ) T describe, respectively, the output and input waves. The off-diagonal elements of S correspond to the transmission coefficients from left and right which, here, are equal due to reciprocity (S = S T is symmetric with respect to the diagonal). The reflection coefficients from ports 1 and 2 are noted respectively r 1 and r 2 .

In this section, we consider networks of identical onedimensional waveguide segments with the same length L, such as the ones shown in Fig. 1(b-d). We note that networks of one-dimensional waveguides are known as "quantum graphs" and have been studied extensively in the last decades; an excellent introduction to the field is given in [START_REF] Berkolaiko | Introduction to Quantum Graphs[END_REF]. Experimentally, they could be realized, for instance, in the form of microwave networks [START_REF] Hul | Experimental simulation of quantum graphs by microwave networks[END_REF][START_REF] Hofmann | Spectral duality in graphs and microwave networks[END_REF] or networks of thin acoustic waveguides [START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF][START_REF] Coutant | Topological two-dimensional Su-Schrieffer-Heeger analog acoustic networks: Total reflection at corners and corner induced modes[END_REF][START_REF] Coutant | Acoustic Su-Schrieffer-Heeger lattice: Direct mapping of acoustic waveguides to the Su-Schrieffer-Heeger model[END_REF]; in this paper, we will focus on the former for our experimental validation, while we focus on the latter for our numerical results.

Let us now investigate the scattering off such a waveguide network, and focus on the behavior of the reflection b-d): Different waveguide networks (upper panel) with their respective frequency-dependent reflection coefficients r1(f ) (solid line), r2(f ) (dashed line) in the complex plane shown in the lower subpanel. These reflection coefficients have been computed using the method from [START_REF] Kottos | Quantum graphs: A simple model for chaotic scattering[END_REF] for acoustic waves with visco-thermal losses (Eqs. (E1) and (E2) with loss-coefficient α = 3•10 -5 ), with a frequency f between 0 and 1000 Hz. Each waveguide has a length L = 0.1m and a diameter R = 2mm. coefficients r 1 (from left) and r 2 (from right). To do so, here, we consider 1D acoustic waves of frequency ω = c 0 k where c 0 is the sound speed and k the wavenumber. The reflection of the two-port system is given by the scattering matrix calculated using the 1D Helmholtz equation in each waveguide and the conservation of acoustic flux at the connections [START_REF] Kottos | Quantum graphs: A simple model for chaotic scattering[END_REF][START_REF] Richoux | Multi-functional resonant acoustic wave router[END_REF]; an explicit expression for the scattering matrix is given in Appendix A. Note that this 1D approximation is valid for acoustic propagation through the waveguides, assuming that L ≫ w, where L (w) is the length (width) of each waveguide.

For a mirror-symmetric network, we expect equireflectionality, that is, r 1 = r 2 , simply due to the total symmetry of the scattering problem. For illustration purposes, we show the two reflection coefficients of the network of Fig. 1(b) in the complex plane in the bottom panel, verifying equireflectionality. Notice that here we have considered the effect of homogeneously distributed losses since the reflection coefficients are inside the unit circle (see below for details). Evidently, such losses maintain the mirror symmetry.

Breaking this mirror symmetry, as is the case in Fig. 1(c), is expected to destroy equireflectionality. That is, r 1 ̸ = r 2 except for some special frequencies where the two reflection coefficients coincide. This is clearly seen in the bottom panel in Fig. 1(c). On the other hand, for a special asymmetric setup as the one in Fig. 1(d) we find that the system is equireflectional, i.e., r 1 = r 2 , as shown in the bottom panel, even in the presence of losses. This surprising result is not obtained by chance or coincidence. In fact, below, we show how such asymmetric equireflectional networks can be designed, and how their equireflectionality can be explained through the recently introduced concept of latent symmetry. The acoustic pressure field here corresponds to a symmetric input (same amplitude and phase) from both ports at f = 625Hz.

B. Review of latent symmetry in waveguide networks

To understand the highly symmetric scattering properties of the setup depicted in Fig. 1(d), we first need to study the eigenmodes of the corresponding cavity that is obtained by closing the connections to the leads, as shown in Fig. 2.

Before we start our theoretical treatment, let us briefly inspect Fig. 2 (a) in more detail. There, we show the wavefield of one of the cavity's eigenmodes. As can be seen, this mode is, to a good approximation, a planewave mode: While the wave field changes along a given waveguide (longitudinal direction), it is constant along its transverse direction. This is, of course, expected from the fact that our narrow waveguides are effectively onedimensional only.

What is interesting about the eigenmode ϕ depicted in Fig. 2 (a) is that its amplitudes at the two points a and b are the same. This is no coincidence. Indeed, as we shall show below, all eigenmodes of the cavity fulfill either ϕ(a) = ϕ(b) or ϕ(a) = -ϕ(b). In other words, the eigenmodes show point-wise parity at the two points a, b. We note that this property for acoustic waveguides has been studied in detail in Ref. [START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF] and here we review the basic results.

The basic idea for finding the eigenmodes of the waveguide network is then as follows: Since our waveguides are effectively one-dimensional, knowledge of the wave field only at at the junctions and at end point (like a, b) of the network is sufficient to reconstruct the wavefield in the entire network. Thus, we can simplify the problem by discretizing the network, as shown in the inset of Fig. 2 (a). The discrete model is nothing but a graph: A collection of points (vertices)-each of them corresponding to a junction or to a waveguide endpoint-interconnected by lines (edges), each of which representing a waveguide.

Following the above idea of discretizing the problem, as sketched above, it can be shown that the problem of finding the eigenmodes ϕ of the network is equivalent to the following generalized eigenvalue problem

A ϕ = cos(kL)B ϕ (2)
with cos(kL) being the eigenvalue, the N -dimensional eigenvector ϕ denoting the values of the eigenmode ϕ at the N vertices. The matrix A describes the topology of the setup, with A i,j = 1 if the vertices i, j are connected by a waveguide, and A i,j = 0 otherwise. The matrix B is diagonal, with B i,i = j A i,j . We then introduce the "Hamiltonian"

H = √ B -1 A √ B -1 , which
arises naturally from Eq. ( 2) through the transformation y = √ Bϕ. We are then left with the eigenvalue problem Hy = cos(kL)y .

(

) 3
Solving Eq. ( 3) gives us the natural frequencies of the eigenmodes for the network cavity, as well as the values of the field at the nodes.

The claim is that, there are two different ways [START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF] to obtain point-wise parity at some points a, b for all eigenmodes. The first is through geometrical symmetry, that is, when the system is invariant under a symmetry operation-usually a reflection, as is the case for Fig. 1 (a)-which maps the junctions a and b onto each other. The second and non-obvious case occurs when the setup is not invariant under a geometrical operation that maps a and b onto each other. This is the case for the network of Fig. 2 (a) and Fig. 1 (d). In fact, the corresponding discrete network, as shown in the inset of Fig. 2 (a), has been designed such that the following relation of the matrix powers of H is fulfilled

(H m ) a,a = (H m ) b,b (4)
for m = 1, . . . , N . This has been found to ensure the pointwise parity of a and b for all the modes [START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF]. Note that Eq. ( 4) is automatically satisfied for a mirror symmetric setup. As we demonstrate in Appendix D, the relations Eq. ( 4) can also be translated into a set of geometric rules that can be visually checked. The simplest of these rules says that the number of neighbors of points a and b-that is, the number of waveguide that are directly connected to these points-have to be the same.

For the network of Fig. 2 (a), both a and b have only one neighbor.

Interestingly, the relations Eq. ( 4) have also a consequence on eigenvalues: If a matrix H satifies these, then the two matrices H \ a and H \ b-obtained from H by removing the ath (or bth) row and column-have the same eigenvalue spectra. We thus say that a and b are cospectral.

Before continuing, let us mention that the relations Eq. ( 4) lead to another very interesting consequence. By performing a certain dimensional reduction onto the two sites {a, b}-the so-called isospectral reduction [START_REF] Bunimovich | Isospectral Transformations: A New Approach to Analyzing Multidimensional Systems and Networks[END_REF], the resulting reduced system can be described by an effective Hamiltonian with a reflection symmetry. This explains the point-wise parity of eigenvectors on a, b. Since this symmetry becomes in general only apparent after the dimensional reduction, we call the Hamiltonian H latently reflection symmetric [START_REF] Smith | Hidden symmetries in real and theoretical networks[END_REF][START_REF] Röntgen | Latent symmetry induced degeneracies[END_REF][START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF].

C. Equireflectionality in Acoustic Waveguide Networks

Let us now investigate the scattering properties of the system by connecting the points a and b to two waveguides with the same cross-section, as shown in Fig. 2(b). For this setup, the two incoming plane waves can be related to the outgoing waves through Eq. (1). To find the corresponding scattering coefficients we first relate the acoustic pressure and its derivative between points a and b through a matrix involving the eigenmodes of the closed cavity. To do so, we use the integral formalism with an expression of the Green function of the closed cavity. Then, as explained in detail in section IV of the supplemental material of [START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF], it can be shown that the former matrix is mirror symmetric due to the point-wise parity of the eigenmodes. This leads to equireflerectionality r 1 = r 2 of the corresponding scattering problem and thus the surprising result of Fig. 1 (d).

So far, we have focussed on the case of thin waveguides. When deviating from this limiting case of w ≪ L, 3D-effects at the junctions of the waveguides will play a role due to evanescent waves; their severity is, clearly, system-dependent (number of junctions, geometry of the connection etc). We calculated the reflection coefficients for various values of w/L and show some characteristic results in Fig. 3. In particular we show the reflections r 1 (f ), r 2 (f ) for the latently-symmetric scattering network of Fig. 1 (d). As can be seen, for this setup, r 1 = r 2 holds quite well even when w/L is far from the limiting value zero, although some deviations are seen in Fig. 3 (b). In Fig. 3(c) we show a more quantitative comparison of the behavior of these two reflections by showing the absolute value of the frequency-dependent difference r 1 (f ) -r 2 (f ) for different sidelengths w. In all of the calculations for Fig. 3, we have used viscothermal losses, which, as can be seen, do not alter equireflectionality. Indeed, it turns out that this kind of losses does not at all impact the latent symmetry of the cavity or the equireflectionality of the (open) system (see Appendix C for details).

We stress that the network depicted in Fig. 1 (d) is, by far, not the only one with latent symmetry-induced equireflectionality. Indeed, as we show in the following, one can easily generate a plethora of such networks.

Perhaps the easiest approach is to generate a large number of different networks and test each of them for latent symmetry. Opening the network at two latently symmetric points will then lead to equireflectional scattering.

Crucially, both of these tasks-generating networks and testing for latent symmetry-can be efficiently done. For the first task, it is important to remember the following: Since we work with strictly one-dimensional (singlemode) waveguides, the actual geometry of the network is not of importance. Instead, what matters is only how the waveguides are interconnected; that is, the network's topology. Mathematically speaking, our waveguide networks are nothing but so-called simple graphs, which are extensively treated in graph theory [START_REF] Newman | Networks: An Introduction[END_REF]. This conceptual overlap allows us to lend tools from graph theory; in particular, we can rely on the powerful nauty suite [START_REF] Mckay | Practical graph isomorphism, II[END_REF]. With nauty, generating a large number of networks is rather efficient; for instance, generating all 11,716,571 distinct waveguide networks with N = 10 sites takes less than 10 seconds on a modern laptop [30]. For the second step, we need to check each of the obtained networks for latent symmetries. To this end, one needs to build the matrices A, B, H. Checking whether the network has a latent symmetry then boils down to computing the first N -1 matrix powers and checking whether there are sites a, b that fulfill Eq. ( 4). In this manner, a plethora of possible networks can be found. Some examples are shown in Fig. 4.

D. Experimental Observation in Microwave

Transmission-Line Networks

As mentioned earlier, besides acoustic waveguide networks, microwave transmission-line networks lend themselves well to experimentally observe the discussed scattering signature of latent symmetry, namely equireflec- tion. Indeed, as can be easily shown [START_REF] Hul | Experimental simulation of quantum graphs by microwave networks[END_REF][START_REF] Hofmann | Spectral duality in graphs and microwave networks[END_REF][START_REF] Mubarack | Transmission Lines, Quantum Graphs and Fluctuations on Complex Networks[END_REF], the mathematical considerations in Section II B can be directly applied to transmission-line networks, if only we replace the pressure field p with the voltage U . The only difference between the two platforms is that, in acoustics, a waveguide with a closed end (acoustic hard wall) features Neumann boundary conditions on that end; to achieve the same in transmission line networks, we need to make the corresponding cable open-ended. In summary, and taking into account this correspondence, the structures with latent symmetry from Fig. 1(d) and Fig. 4 are directly applicable to a microwave realization.

We have hence built the networks from Fig. 1(d) and Fig. 4(a) using 50 cm-long coaxial cables and measured their scattering parameters with a vector network analyzer (VNA, Rhode & Schwarz ZVA 67, 10 MHz -67 GHz). The finite propagation delay in the junctions is equivalent to the case of point-like junctions with slightly longer waveguides. The latter is in line with our theoretical model. We also make sure that the effective length of each waveguide is identical.

We observe in Fig. 5 excellent agreement of the complex-valued reflection coefficients (i.e., in terms of both magnitude and phase) with the expected equireflection condition at all frequencies. Moreover, a clear periodicity of the reflection spectrum is visible. Furthermore, it is apparent that the absorption strength in-Figure 4. Different waveguide networks whose scattering matrix is equireflectional. Just as in Fig. 1, the reflection coefficents are computed for each network and plotted (in the complex plane) below each network. Since the networks are equireflectional, the two curves r1(f ) (solid line) and r2(f ) (dashed line) are identical. The reflection coefficents have been computed with the same parameters as in Fig. 1, though now up to f = 2000 Hz.

creases monotonously with frequency, as expected.

III. POWERS OF THE SCATTERING MATRIX AND PERFECT ABSORPTION

Although two-port systems provide an intuitive and simple setup to study wave scattering, many applications are actually based on N -port configurations (circulators, splitters, multiplexers etc.). In view of these implementations, in this Section, we study latent symmetry and generalisations of it directly on the scattering matrix itself and not on the underlying structure. In such a way the results obtained here do not depend on the specific physical system and are rather general for any N-port scattering matrix S. We first note that any reciprocal scattering matrix is symmetric S = S T . Let us then consider a case where this matrix is also latently symmetry; that is, it fulfills

(S m ) a,a = (S m ) b,b , ∀ m = 1, . . . , N (5)
for some a and b. Note that this relation for m = 1 implies that ports a and b have the same reflection coefficient. We emphasize that the connection between matrix powers, cospectrality and eigenvectors (presented in the last section) is valid for the scattering matrix as well. Thus, if Eq. ( 5) is valid, the eigenvectors of the S matrix will have parity between the elements a and b. However, the eigenvectors of the scattering matrix, except for particular cases (as we show below) are not so useful. We remind the reader that Eq. ( 5) also implies cospectrality between S \ a and S \ b but gives no information about the eigenvalues of S itself. To bypass this, let us go one step further and generalise Eq. ( 5) by introducing the following relation

(S m ) a,a = c • (S m ) b,b ∀ m = 1, . . . , N . (6)
We call this relation, Eq. ( 6), scaled cospectrality of a matrix S with a scaling factor c. Interestingly, this rather simple generalisation has direct impact on the eigenvalues of the matrix S: whenever c ̸ = 1, scaled cospectrality enforces the existence of one zero eigenvalue (more zeros might appear accidentally). To show why this is the case, let us start by realizing that, from the Cayley-Hamilton theorem, we know that S N can be written as a polynomial in the first N -1 powers of S, that is,

S N = - N -1 m=0 a m S m , (7)
where S 0 ≡ I is the identity matrix and where the coefficients a m are taken from the characteristic polynomial of S, given by P (x) = det (x I -S) = N m=0 a m x m . Now, since the matrix S fulfills S N a,a = c S N b,b , it follows from Eq. ( 7) that

-a 0 I a,a - N -1 m=1 a m (S m ) a,a = c -a 0 I b,b - N -1 m=1 a m (S m ) b,b (8)
It follows that a 0 = 0 whenever c ̸ = 1 and since a 0 is proportional to the determinant of S, we see that S has to have a vanishing eigenvalue. The existence of a zero eigenvalue of S directly implies that there must be a monochromatic adapted wavefront that can be injected at the real frequency at which S has the zero eigenvalue such that it is perfectly absorbed within the system [START_REF] Chen | Perfect absorption in complex scattering systems with or without hidden symmetries[END_REF]. This CPA requires the system to have a finite amount of absorption loss and hence a sub-unitary S matrix. Indeed, in the absence of any absorption loss, the S matrix would be unitary and the magnitude of its determinant would be unity, implying that a 0 cannot be zero because a 0 is proportional to det(S). Let us remark for the interested reader that in the literature there is also another generalization of cospectrality, the so-called fractional cospectrality [START_REF] Chan | Fundamentals of fractional revival in graphs[END_REF], which, however, does not enforce the presence of a zero eigenvalue.

In this section, we have so far not explicitly stated the frequency dependence of S. However, a scattering system can in general only fulfill the conditions for scaled cospectrality at discrete frequencies, which is hence in line with the fact that the zeros of S only occur at discrete (potentially complex) frequencies. In other words, if one optimizes the scattering system such that it has scaled cospectrality (and thus CPA) at some frequency f 0 , there is no reason to expect scaled cospectrality at any frequency other than f 0 .

In the remainder of this section, we will show that systems with scaled cospectrality allow us to achieve CPA with wavefronts that have a customized imbalance in terms of the weights of different channels. Before explaining the origin of this feature in detail, let us briefly contextualize this feature within the recent literature on achieving CPA. The S matrix of a generic (arbitrarily complex) scattering system does not necessarily have a zero eigenvalue, but can be tuned to have one [START_REF] Pichler | Random antilasing through coherent perfect absorption in a disordered medium[END_REF][START_REF] Chen | Perfect absorption in complex scattering systems with or without hidden symmetries[END_REF][START_REF] Fyodorov | Distribution of zeros of the S-matrix of chaotic cavities with localized losses and coherent perfect absorption: Non-perturbative results[END_REF][START_REF] Li | Random matrix theory approach to chaotic coherent perfect absorbers[END_REF]. With sufficient tunable degrees of freedom, the real frequency at which the zero eigenvalue occurs can be controlled, too [START_REF] Imani | Perfect absorption in a disordered medium with programmable meta-atom inclusions[END_REF][START_REF] Frazier | Wavefront shaping with a tunable metasurface: Creating cold spots and coherent perfect absorption at arbitrary frequencies[END_REF][START_REF] Del Hougne | On-Demand Coherent Perfect Absorption in Complex Scattering Systems: Time Delay Divergence and Enhanced Sensitivity to Perturbations[END_REF]. However, in structures without symmetry the required CPA wavefront is in general highly asymmetrical and difficult to generate [38]. This makes it more challenging to realize the envisaged use of CPA to interferometrically control light with light, and hence without any non-linearity [START_REF] Zhang | Controlling light-with-light without nonlinearity[END_REF][START_REF] Baranov | Coherent perfect absorbers: linear control of light with light[END_REF]. A generic way to impose CPA at a desired frequency and with an arbitary CPA wavefront was demonstrated in Ref. [START_REF] Del Hougne | Coherent wave control in complex media with arbitrary wavefronts[END_REF] by tuning a massively parametrized chaotic cavity. Here, we provide an alternative route through scaled cospectrality to obtain structures featuring CPA whose wavefront has a prescribed imbalance. Thereby, only a weak control wave can modulate a strong signal, provided that the necessary phase and amplitude relation between the two is respected.

We now illustrate the above results with an implementation of the scaled cospectrality using a particular example of a 3-port acoustic network with lossy acoustic waveguides (see Appendix C for details). For simplicity we choose to work with an extension of the 2-port system of the previous section by symmetrically coupling an arbitrary structure to a and b as shown in Fig. 6 (a). The symmetric connection and the added structure are highlighted by light gray and dark gray respectively.

The two latently symmetric points are symmetrically connected to ports 1 and 2, and thus the scattering matrix of the 3-port network reads

S =   r t t ′ t r t ′ t ′ t ′ r ′   . (10)
with r, r ′ denoting reflection coefficients, and with t, t ′ transmission coefficients. Evidently the proposed design ensures the equireflectionality (r) of the two symmetric ports 1 and 2 while an additional reflection (r ′ ) from port 3 is introduced. Due to reciprocity and symmetry we are thus left with only two different transmission coefficients t and t ′ . For this particular matrix we impose scaled cospectrality between ports 1 (or equivalently 2) and 3 by demanding that it satisfies Eq. ( 6). Note that for m = 1 we simply get a scaling between the reflections, r = cr ′ . Then, solving for the higher powers we obtain two possible solutions which are discussed in the following.

CPA without any zero-input channel: One solution requires t = -r ′ /2, and To illustrate the effects of the zero eigenvalue of S and its relation to CPA, we define the following eigenvalue problem

t ′ = r ′ √ 2c -1/2 . ( 11
Sx = λx. (12)
When S satisfies Eq. ( 11), the corresponding eigenvalues of the scattering matrix are λ 1 = (0, (c + 1/2)r ′ , (c + 1/2)r ′ ). As expected by our definition of the generalised cospectrality, the scattering matrix acquires a zero eigenvalue. This zero eigenvalue implies that if one uses as an input the corresponding eigenvector, the output is zero and the input wavefront is perfectly absorbed. In fact, here the corresponding eigenvector is

x 1 = (-1/ √ 2c -1, -1/ √ 2c -1, 1) T . (13)
Importantly, the imbalance of the input wavefront is solely controlled by the scaling factor c which we can prescribe.

In order to validate our results, we once again consider the case of airborne acoustic waveguides including viscothermal losses. The elements of the S matrix in such a case can be calculated using different techniques such as the star product (see Appendix E for details). We optimize the lengths of several parts of the network (while maintaining the symmetry between ports 1 and 2) such that Eq. ( 11) is satisfied for a predefined value of c and a suitable frequency f 1 . To illustrate our result, we assume an input to the 3-port network in the form of a vector p in = x 1 and thus the output vector is naturally given by p out = S p in . Using the output we can calculate absorption from such networks as

A ≡ 1 - |p out | 2 |p in | 2 . (14)
The absorption for two different networks is shown in Fig. 6(b). The black thick (orange thin) line corresponds to a network satisfying Eq. ( 11) with c = 100 (c = 0.01).

The details of the geometry of the obtained networks are given in Appendix E. According to our design, at the prescribed frequency f 1 , CPA of the input wave is achieved, indicated by A = 1. We additionally observe that with relatively small variations of the lengths of the waveguides (see Appendix E), we achieve CPA for a very different input vector of an amplitude 10 times larger (smaller) for port 1/2 compared to port 3. CPA with one zero-input channel: The other possible solution for S to acquire scaled cospectrality is when

t = r = cr ′ and t ′ = √ cr ′ , (15)
with the eigenvalues of S then being λ 2 = (0, 0, (2c+1)r ′ ). The corresponding eigenvectors of the two zero eigenvalues are x

(1) 2 = (1, -1, 0) T and x (2) 2 = (1, 1, -2 √ c)
T . This case is even more interesting since it features a two-fold degenerate zero eigenvalue [START_REF] Piper | Total absorption by de-generate critical coupling[END_REF]. Consequently, for any choice of α, β, inputs of the form αx will be perfectly absorbed. Among all possible inputs, one can find the following highly asymmetric one:

x 2 = (-1, 0, √ c) T . (16)
Thus, the 3-port network is able to completely absorb waves with non-zero inputs on only two of its ports and with a relative input between the two active ports prescribed by the scale factor c. We have constructed two such networks using c = 100 and c = 0.01, satisfying Eq. ( 15) at a prescribed frequency f 2 and the corresponding absorption is shown in Fig. 6(c) using the input p in = x 2 . The networks feature CPA (A = 1) at the desired frequency f 2 with non-zero inputs in only two of the three channels. On top of that, the imbalance of the CPA wavefront introduced by the scale factor c results in an almost single sided (or one port) CPA.

IV. CONCLUSIONS

We have studied the scattering properties of various asymmetric waveguide networks and demonstrated that a certain family of networks possesses the same scattering properties as mirror-symmetric ones, i.e. they show equireflectionality. This counterintuitive property of the asymmetric networks stems from a hidden mirror symmetry called latent symmetry. We have validated this finding numerically for acoustic waveguide networks and experimentally for microwave transmission-line networks.

While latent symmetry is mathematically equivalent to certain relations obeyed by the matrix powers, here we generalized these relations and applied them to a generic scattering matrix of an N -port system. The new relations, named scaled cospectrality, were then used to construct networks featuring CPA. Specifically, a scaling factor was used to design systems able to completely absorb wavefronts with prescribed imbalance, a capability that can enable the control of light with a very weak coherent control signal. Overall, our work demonstrates that scattering problems may greatly profit from the notion of latent symmetry-which has very recently been experimentally tested also in non-reciprocal networks [START_REF] Sol | Covert Scattering Control in Metamaterials with Non-Locally Encoded Hidden Symmetry[END_REF]-and of matrix power relations in general.

with

k = 2πf c s + (1 + i) α √ f R ; (E2)
where α is the loss coefficient (see below) [START_REF] Kosten | Sound Absorbing Materials[END_REF], R = 1 cm the waveguide diameter, with c s = 343 m/s the velocity of sound and f being the frequency. Using the continuity of the pressure and the conservation of the flux at each junction, we then obtain the scattering matrix S, depending only on the frequency f and the lengths L m . The optimized lengths (in meters, rounded to 4 digits after the decimal point) are shown in the table below. The cases a and b correspond to the scattering matrix fulfilling Eq. [START_REF] Röntgen | Latent symmetry induced degeneracies[END_REF] and Eq. ( 15) of the main text, respectively. case a, c = 1/100 case a, c = 100 case b, c = 1/100 case b, c = 100

L

Figure 1 .

 1 Figure 1. (a) A generic two-port scattering setup. (b-d): Different waveguide networks (upper panel) with their respective frequency-dependent reflection coefficients r1(f ) (solid line), r2(f ) (dashed line) in the complex plane shown in the lower subpanel. These reflection coefficients have been computed using the method from [21] for acoustic waves with visco-thermal losses (Eqs. (E1) and (E2) with loss-coefficient α = 3•10 -5 ), with a frequency f between 0 and 1000 Hz. Each waveguide has a length L = 0.1m and a diameter R = 2mm.

Figure 2 .

 2 Figure 2. The network of Fig. 1 (d), realized with thin, square-shaped acoustic waveguides of length L = 0.1m and side length w = 10mm (see Section II C for details). (a) The 10th mode of the acoustic network with closed ends. The inset shows the corresponding discrete model (graph). (b) Schematic of the scattering process by the 2-port network.The acoustic pressure field here corresponds to a symmetric input (same amplitude and phase) from both ports at f = 625Hz.

Figure 3 .

 3 Figure 3. Results from 3D finite-element simulations of scattering off the structure of Fig. 1(d) using COMSOL, with visco-thermal losses modelled through the "narrow region acoustics" function. The waveguides are identical with a length of L = 0.1m, and with a quadratic cross-section with side length w. (a) and (b) show a comparison of r1 (solid line) and r2 (dashed line) for waveguide side lengths w = 5mm and w = 15mm, respectively. The frequency range is between 0 and 2000Hz, which comprises a little more than the first 12 eigenmodes of the underlying cavity. (c) Comparison of |r1 -r2| for three values of w.

Figure 5 .

 5 Figure 5. Experimental scattering measurements on microwave transmission-line networks with latent symmetry. (a) Photography of the experimental setup for the network from Fig. 1(d). (b) and (c) Magnitude and phase of the measured reflection coefficients for the networks from Fig. 1(d) and Fig. 4(a), respectively.

  from which we get a 0 (I) a,a = c a 0 (I) b,b .

Figure 6 .

 6 Figure 6. (a) A symmetric extension of the 2-port setup shown Fig. 1 (d). For (b) and (c), we finetune the lengths of the acoustic waveguides of the network (a) to fulfill scaled cospectrality between ports 1 and 3 with the scattering matrix fulfilling either Eq. (11) [for (b)] or Eq. (15) [for (c)]. For each (b) and (c), we show the frequency-dependent absorption A of the network for the two different choices c = 100 and c = 1/100.

Figure 7 .

 7 Figure 7. Reference figure showing the assignment of the waveguide lengths L1, . . . , L5.

ACKNOWLEDGMENTS

The authors are thankful to M. Pyzh and V. Pagneux for valuable discussions. V.A. acknowledges financial support from the NoHeNA project funded under the program Etoiles Montantes of the Region Pays de la Loire. V. A. is supported by the EU H2020 ERC StG "NASA" Grant Agreement No. 101077954.

Appendix A: Formula for the scattering matrix

In the following, we describe how the scattering matrix of the waveguide networks (assuming one-dimensional waveguides) can be computed. The expressions are taken from [START_REF] Kottos | Quantum graphs: A simple model for chaotic scattering[END_REF]. Technically, they rely on describing the network in terms of a so-called quantum graph with Neumann vertex conditions.

Given the topology matrix A describing the waveguide network (see main text), and the diagonal matrix B with entries B i,i = j A i,j , we define the two matrices

and the 2 × N coupling matrix W whose elements are all zero, except W 1,a = W 2,b = 1. L denotes the length of the individual waveguides, and k is the (complex) wavenumber.

The scattering matrix can be shown to be given by [START_REF] Kottos | Quantum graphs: A simple model for chaotic scattering[END_REF]

where I is the 2 × 2 identity matrix. (B1)

We assume that S has no degenerate eigenvalues, which in particular implies that S is diagonalizable. Then, since S = S T , one can normalize its eigenvectors |ϕ i ⟩ such that ⟨ϕ * j |ϕ i ⟩ = δ i,j with the star denoting the complex conjugate. As a consequence, we have 1 = i |ϕ i ⟩ ⟨ϕ * i |. Equipped with this identity, we get

with λ i being the eigenvalue of |ϕ i ⟩. Since the complex exponentials are linearly independent, one can evaluate the above equation independently for each distinct λ i , so that it automatically follows that

Since ⟨ϕ * j |a⟩ = ⟨a|ϕ j ⟩ for all eigenvectors, we obtain

Appendix C: Equireflectionality in a lossy system

In the main text of this work, we considered the case of latent reflection symmetry in a system without losses, and we mentioned that the point-wise parity of eigenmodes remains valid for the case of losses in the one-dimensional waveguides described by a complex velocity c = c r + ic i . In the following, we give the justification for this statement.

To this end, let us assume that we have a network which, in the absence of losses, has a latent reflection symmetry between two junctions a, b. As stated in the main text, for low-enough frequencies, the eigenmodes of this network can be found from the generalized eigenvalue problem

where ϕ contains the pressure of the eigenmode ϕ at the N junctions of the network, and with k = ω/c = ω/c r ≡ k r . The eigenvalues cos(kL) of this problem are completely real and due to latent symmetry, all eigenmodes fulfill p a = ±p b .

Let us then take the usual route for introducing losses: That is, we start from the lossless case-as described through Eq. (C1)-and then let the velocity c = c r + ic i (due to, in particular, thermo-viscous boundary at the surface of the one-dimensional waveguides [START_REF] Kosten | Sound Absorbing Materials[END_REF]). Obviously, this changes the relation between k and ω, but it does not change the eigenvalues cos(kL) or the eigenvectors ϕ. In other words, each eigenvector ϕ with eigenvalue cos(kL) of the lossless system will still be an eigenvector of the lossy system, with unchanged eigenvalue cos(kL). In particular, point-wise parity of eigenmodes is preserved. What changes is the frequency ω corresponding to cos(kL): Since c is complex while k is real (as imposed from the fact that cos(kL) is real), ω is in general complex as well, with the imaginary part being related to the lifetime of this lossy eigenmode.

Appendix D: Using the matrix power relations for deriving a better understanding of latently symmetric waveguide networks

The aim of this section is to showcase a set of intuitive and easily interpretable equivalent conditions that the relations Eq. ( 4)-latent symmetry, that is-impose on a waveguide network. These conditions can be derived by analyzing Eq. ( 4) order by order and subsequently using the relation H = B -1/2 AB -1/2 to derive conditions on the matrices A, B which are directly describing the underlying waveguide network. We restrict ourselves to the first few orders of Eq. ( 4), for which the corresponding equivalent conditions have been derived in [START_REF] Röntgen | Hidden Symmetries in Acoustic Wave Systems[END_REF]; we repeat them here for self-containedness of the present manuscript.

The conditions are then as follows. Firstly, a and b have to have the same number of neighbors. Moreover, in the special case where the number of next-neighbors of a, b is equal to unity, (i) Eq. ( 4) holds for m = 2 if and only if a, b have the same number of next-neighbors, and (ii) Eq. ( 4) holds for m = 3 if and only if

In this relation, N 2 (i) denotes the set of next-neighbors of i, and |N (i)| denotes the degree of site i, that is, the number of neighbors of i. Equipped with the above, we can now analyse the difference between the setups shown in Fig. 1 (c) and Fig. 1 (d) in more detail. To this end, let us close these two systems on the entry-points of the two ports, and call these points a (on the left-hand side of the setup) and b (on the right-hand ride of the setup); compare Fig. 2.

For the asymmetric network of Fig. 1 (c), the two next-neighbors of a have degrees 3 and 2, while the two nextneighbors of b both have degree 3. Thus, the equation Eq. (D1) is not fulfilled and, as a consequence, the relations Eq. ( 4) are not fulfilled for m = 3. On the other hand, the setup of Fig. 1 (d) features an additional waveguide on the left, which equalized the two sides of Eq. (D1). In order to design a system whose scattering matrix has the form of Eq. ( 10) and which fulfills either Eq. ( 15) or Eq. ( 11), we proceed as follows: We started with the system depicted in Fig. 6 (a), whose scattering matrix has the structure of Eq. [START_REF] Bunimovich | Isospectral Transformations: A New Approach to Analyzing Multidimensional Systems and Networks[END_REF]. We then optimized the lengths L 1 , . . . , L 5 (see Fig. 7) and the frequency f in the range between 1000 and 2000 Hz such that Eqs. ( 11) and ( 15) are fulfilled, respectively, at a frequency f . Technically, we modeled the system as a network of ideal waveguides of length L m , each with a transmission coefficient

Appendix F: Confirming the reciprocity of the transmission line networks

In the main text, we have already shown the reflection coefficients for the transmission line measurements of the networks shown in Fig. 1(d) and Fig. 4(a). In Fig. 8, we also show the transmission coefficients, that is, the matrix elements S 1,2 and S 2,1 of the scattering matrix S. As can be seen, the networks are reciprocal in their scattering properties, that is, they fulfill S 1,2 = S 2,1 .