I have done a lot more digging on this. The resonant frequencies in a spherical cavity are given by Eq 18 in this paper: https://www.acs.psu.edu/drussell/Publications/Basketball.pdf. They are dependent on zeros of the derivatives of spherical Bessel functions.

I have an example configuration of a cavity of water (c = 1500 m/s) with radius 2.0 cm surrounded by a layer of glass (outer radius 3.2 cm) and then air. Since the smallest zero of the derivatives of spherical Bessel functions is 2.08, the fundamental frequency of this cavity should be 2.08*1500/(2*3.141*0.02) = 24832 Hz ~ 25 kHz. However, simulating this in k-Wave gives a fundamental of ~20 kHz. I have attached the code I am using to simulate this. I have made sure to place the PML outside of the cavity. What am I missing here? Is this a problem regarding reflections from the PML if the absorption is set improperly?

```
% create the computational grid
Nx = 32; % number of grid points in the x (row) direction
Ny = 32; % number of grid points in the y (column) direction
Nz = 32; % number of grid points in the z (frame) direction
dx = 0.2e-2; % grid point spacing in the x direction [m]
dy = 0.2e-2; % grid point spacing in the y direction [m]
dz = 0.2e-2; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
kgrid.dt = 1e-8;
kgrid.Nt = 100000;
% make the spherical water ball surrounded by air
air_speed = 343;
air_density = 1;
interior_speed = 1500;
interior_density = 1000;
exterior_speed = 2500;
exterior_density = 2000;
medium.density = air_density*ones(Nx,Ny,Nz);
medium.sound_speed = air_speed*ones(Nx,Ny,Nz);
interior = makeBall(Nx, Ny, Nz, Nx/2, Ny/2, Nz/2, 10);
exterior = makeBall(Nx, Ny, Nz, Nx/2, Ny/2, Nz/2, 16);
medium.density(logical(exterior)) = exterior_density;
medium.sound_speed(logical(exterior)) = exterior_speed;
medium.density(logical(interior)) = interior_density;
medium.sound_speed(logical(interior)) = interior_speed;
% initial pressure at center of the domain
disc_magnitude = 30; % [Pa]
disc_x_pos = Nx/2; % [grid points]
disc_y_pos = Ny/2; % [grid points]
disc_z_pos = Nz/2; % [grid points]
disc_radius = 5; % [grid points]
disc_2 = disc_magnitude * makeBall(Nx, Ny, Nz, disc_x_pos, disc_y_pos, disc_z_pos, disc_radius);
source.p0 = disc_2;
% sensor mask uniformly distributed over space
sensor.mask = zeros(Nx,Ny,Nz);
sensor.mask(1:5:Nx,1:5:Ny,1:5:Nz)=1;
% run the simulation
input_args = {'DataCast','single','DeviceNum', 2, 'PMLInside', false,'PMLSize',16};
sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor,input_args{:});
% sum over the spacial domain, fft the signal, and plot the frequency domain
FD = abs(fft(sum(sensor_data,1)));
FD = FD/max(FD,[],'all');
freq = (0:(kgrid.Nt-1))*(1/kgrid.dt)/kgrid.Nt;
semilogx(freq(1:end/2),FD(1:end/2));
```

Here is the output I am getting from the above code: https://drive.google.com/file/d/1R5G82cPa7DsZB4attrFqCW7Rqv5hy9Pq/view?usp=sharing

PC