In HFSS, create an eigenmode solution simulation. Create the appropriate boundary conditions (most general, Master/Slave with a 180 (and/or 0) degree phase shift) and solve for the first two modes. This will give you the frequencies of these modes as well as the field distributions.
To add to my previous response, I've created a basic example in HFSS (attached). Please note that this eigenmode solution will NOT take into account the effect of fringing fields, so if you're designing something like a patch antenna the frequency will be slightly off.
I've created a system of two copper conductors sandwiching a slab of 60 mil FR-4. Around the sides of the FR-4, we place the master-slave boundaries, one for each axis. The phase shift entered for PhaseX and PhaseY will determine the modes that are solved for. These boundaries simply match their fields with the specified phase shift. So, for example a PhaseX of 0 degrees and PhaseY if 180 degrees will solve all TM_(even)(odd) modes in the cavity.
The eigenmode simulation is set up so that it will only solve for the first two modes - this is given in order of lowest frequency to highest. So in order to solve for additional higher-order modes, just solve for more modes.
Please ask if you have any further question about this.
The modes in the patch resonator are ideally not affected by the air above the patch (admittedly not very physical, but approximate). Obviously if you are going to simulate something like an antenna, you will need to include the air box.
The eigenmode simulation is for simulating something with mostly hard (i.e., not open) boundary conditions. This translates into a more idealistic setup - for example, my setup won't include the effects of fringing fields, but it will give you a pretty good idea of what to expect in a resonant patch. You could always add in some space between the patch and the boundaries, along with an air box, but don't make it too large as the eigenmode solver might not be able to work with anything that is largely open.
I'm not sure of what you mean by TE-TM plots. Solve for additional modes if you want to find additional TM or TE modes. If you mean plotting more than one mode at a time, you can do this by setting the excitation strength of each mode under "Field Overlays" (In the left-hand menu) -> "Edit Sources".