## 3 Materials and Methods  
### Simulations  
The k-Wave Matlab toolbox was used to simulate 2D non-linear wave propagation through an inhomogeneous dispersive medium [19; 20]. A CoMTUS system formed by two identical linear arrays, similar to the ones experimentally available, was simulated as follows. Each of the arrays had a central frequency of 3 MHz and 144 active elements in both transmit and receive, with an element pitch of 240 \(\upmu\)m and kerf of 40 \(\upmu\)m. For PWs, the modelled transducer had an axial focus of infinity with all 144 elements firing simultaneously. A simulation was performed for each transmit event, i.e., each PW at a certain angle. In total, seven transmit simulations per linear array were performed to produce a PW dataset, which covered a total sector angle of 30\({}^{\circ}\) (from \(-15^{\circ}\) to \(15^{\circ},5^{\circ}\) step). In the case of CoMTUS, this results in 14 transmit events in total (seven PWs per array). This PW sequence was chosen to match in resolution a focused system with an F-number of 1.9, decimating the required number of angles by a factor of six to optimize the simulation time without affecting resolution [21]. The computational grid was set following convergence testing from a previous study [22]. The spatial grid was fixed at 40 \(\upmu\)m (six grid points per wavelength) with a time step of 1.3 ns. The power law exponent to model frequency-dependent attenuation was set to 1.5 [23]. Received signals were downsampled at 30.8 MHz. Channel noise was introduced to the RF simulated data as white Gaussian noise. The signal-to-noise ratio (SNR) for all simulated data was set to 35 dB at a 50 mm imaging depth. The US pulses were propagated through heterogeneous scattering media using tissue maps defined by the SOS, density, attenuation and nonlinearity (see Table 1). A medium defined only with the properties of general soft tissue was used as the control case. To model the scattering properties observed in vivo, scatterers with a size smaller than the wavelength of the transmitted signal were added to the tissue maps. In order to fully develop speckle [24], a total of 15 scatterers of 40 \(\upmu\)m in diameter, with random spatial position and amplitude (defined by a 5% difference in SOS and density from the surrounding medium), were added per resolution cell, the size of which was assumed to be 0.5 mm\({}^{2}\). Three point-like targets and an anechoic lesion were included in the media to allow the measurement of the basis metrics for comparing the imaging quality for different scenarios. A circular anechoic lesion of 12 mm in diameter located at the centre of the aperture of the common FOV was modelled as a region without scatterers. The point-like targets were simulated as circles of 0.2 mm in diameter with a 25% difference in SOS and density with the surrounding tissue to generate appreciable reflection. The same realization of scatterers was superimposed on all maps and through the different simulations to keep the speckle pattern in the CoMTUS system, so any changes in the quality imaging metrics were due to changes in the overlying tissues, the imaging depth and the acoustical field. The k-Wave Matlab toolbox uses a Fourier co-location method to compute spatial derivatives and numerically solve the governing model equations, which requires discretisation of the simulation domain into an orthogonal grid. Consequently, continuously defined acoustic sources and media need to be sampled on this computational grid, introducing staircasing errors when sources do not exactly align with the simulation grid [25; 26]. To minimize these staircasing errors, the transmit array was always aligned to the computational grid, i.e., simulations were performed in the local coordinate system of the transmit array.