diff --git a/tests/test_IO.py b/tests/test_IO.py
index d880bedca6..c75591c543 100644
--- a/tests/test_IO.py
+++ b/tests/test_IO.py
@@ -45,7 +45,6 @@ def test_simulation_preserve_types():
                 size=(inf, inf, 0),
                 source_time=st,
                 direction="+",
-                polarization="Ex",
                 waist_radius=1,
             ),
         ],
diff --git a/tests/test_components.py b/tests/test_components.py
index 2fe67ff606..58f657c82b 100644
--- a/tests/test_components.py
+++ b/tests/test_components.py
@@ -481,30 +481,26 @@ def test_source_times():
     # c.amp_time(ts)
 
 
-def test_VolumeSource_directional():
+def test_FieldSource():
     g = GaussianPulse(freq0=1, fwidth=0.1)
+    mode = Mode(mode_index=0)
 
     # test we can make planewave
-    s = PlaneWave(size=(0, 1, 1), source_time=g, polarization="Ez", direction="+")
+    s = PlaneWave(size=(0, 1, 1), source_time=g, pol_angle=np.pi / 2, direction="+")
 
-    # test we can make planewave
-    s = GaussianBeam(size=(0, 1, 1), source_time=g, polarization="Ez", direction="+")
+    # test we can make gaussian beam
+    s = GaussianBeam(size=(0, 1, 1), source_time=g, pol_angle=np.pi / 2, direction="+")
 
-    # test that non-planar geometry crashes plane wave
-    with pytest.raises(ValidationError) as e_info:
-        s = PlaneWave(size=(1, 1, 1), source_time=g, polarization="Ez", direction="+")
+    # test we can make mode source
+    s = ModeSource(size=(0, 1, 1), direction="+", source_time=g, mode=mode)
 
     # test that non-planar geometry crashes plane wave and gaussian beam
     with pytest.raises(ValidationError) as e_info:
-        s = PlaneWave(size=(1, 1, 0), source_time=g, polarization="Ez", direction="+")
+        s = PlaneWave(size=(1, 1, 1), source_time=g, pol_angle=np.pi / 2, direction="+")
     with pytest.raises(ValidationError) as e_info:
-        s = GaussianBeam(size=(1, 1, 1), source_time=g, polarization="Ez", direction="+")
-
-
-def test_VolumeSource_modal():
-    g = GaussianPulse(freq0=1, fwidth=0.1)
-    mode = Mode(mode_index=0)
-    m = ModeSource(size=(0, 1, 1), direction="+", source_time=g, mode=mode)
+        s = GaussianBeam(size=(1, 1, 1), source_time=g, pol_angle=np.pi / 2, direction="+")
+    with pytest.raises(ValidationError) as e_info:
+        s = ModeSource(size=(1, 1, 1), direction="+", source_time=g, mode=mode)
 
 
 """ monitors """
diff --git a/tests/test_grid.py b/tests/test_grid.py
index 4165d757c0..1b5a6e7352 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -20,6 +20,7 @@ def test_field_grid():
     c = Coords(x=x, y=y, z=z)
     f = FieldGrid(x=c, y=c, z=c)
 
+
 def test_grid():
 
     boundaries_x = np.arange(-1, 2, 1)
@@ -48,9 +49,9 @@ def test_sim_nonuniform_small():
     grid_size_x = [2, 1, 3]
     sim = td.Simulation(
         center=(1, 0, 0),
-        size=(size_x, 4, 4), 
+        size=(size_x, 4, 4),
         grid_size=(grid_size_x, 1, 1),
-        pml_layers=[td.PML(num_layers=num_layers_pml_x), None, None]
+        pml_layers=[td.PML(num_layers=num_layers_pml_x), None, None],
     )
 
     bound_coords = sim.grid.boundaries.x
@@ -61,13 +62,16 @@ def test_sim_nonuniform_small():
 
     # checks the bounds were adjusted correctly
     # (smaller than sim size as is, but larger than sim size with one dl added on each edge)
-    assert np.sum(dls) <= size_x + num_layers_pml_x*dl_min + num_layers_pml_x*dl_max
-    assert np.sum(dls)+dl_min+dl_max >= size_x + num_layers_pml_x*dl_min + num_layers_pml_x*dl_max
+    assert np.sum(dls) <= size_x + num_layers_pml_x * dl_min + num_layers_pml_x * dl_max
+    assert (
+        np.sum(dls) + dl_min + dl_max
+        >= size_x + num_layers_pml_x * dl_min + num_layers_pml_x * dl_max
+    )
 
     # tests that PMLs were added correctly
     for i in range(num_layers_pml_x):
-        assert np.diff(bound_coords[i:i+2]) == dl_min
-        assert np.diff(bound_coords[-2-i:len(bound_coords)-i]) == dl_max
+        assert np.diff(bound_coords[i : i + 2]) == dl_min
+        assert np.diff(bound_coords[-2 - i : len(bound_coords) - i]) == dl_max
 
     # tests that all the grid sizes are in there
     for size in grid_size_x:
@@ -78,7 +82,8 @@ def test_sim_nonuniform_small():
         assert dl in grid_size_x
 
     # tests that it gives exactly what we expect
-    assert np.all(bound_coords == np.array([-12,-10,-8,-6,-4,-2,0,1,4,7,10,13,16]))
+    assert np.all(bound_coords == np.array([-12, -10, -8, -6, -4, -2, 0, 1, 4, 7, 10, 13, 16]))
+
 
 def test_sim_nonuniform_large():
     # tests when the nonuniform grid extends beyond the simulation size
@@ -88,9 +93,9 @@ def test_sim_nonuniform_large():
     grid_size_x = [2, 3, 4, 1, 2, 1, 3, 1, 2, 3, 4]
     sim = td.Simulation(
         center=(1, 0, 0),
-        size=(size_x, 4, 4), 
+        size=(size_x, 4, 4),
         grid_size=(grid_size_x, 1, 1),
-        pml_layers=[td.PML(num_layers=num_layers_pml_x), None, None]
+        pml_layers=[td.PML(num_layers=num_layers_pml_x), None, None],
     )
 
     bound_coords = sim.grid.boundaries.x
@@ -101,13 +106,16 @@ def test_sim_nonuniform_large():
 
     # checks the bounds were adjusted correctly
     # (smaller than sim size as is, but larger than sim size with one dl added on each edge)
-    assert np.sum(dls) <= size_x + num_layers_pml_x*dl_min + num_layers_pml_x*dl_max
-    assert np.sum(dls)+dl_min+dl_max >= size_x + num_layers_pml_x*dl_min + num_layers_pml_x*dl_max
+    assert np.sum(dls) <= size_x + num_layers_pml_x * dl_min + num_layers_pml_x * dl_max
+    assert (
+        np.sum(dls) + dl_min + dl_max
+        >= size_x + num_layers_pml_x * dl_min + num_layers_pml_x * dl_max
+    )
 
     # tests that PMLs were added correctly
     for i in range(num_layers_pml_x):
-        assert np.diff(bound_coords[i:i+2]) == grid_size_x[0]
-        assert np.diff(bound_coords[-2-i:len(bound_coords)-i]) == grid_size_x[-1]
+        assert np.diff(bound_coords[i : i + 2]) == grid_size_x[0]
+        assert np.diff(bound_coords[-2 - i : len(bound_coords) - i]) == grid_size_x[-1]
 
     # tests that all the grid sizes are in there
     for size in grid_size_x:
diff --git a/tests/test_plugins.py b/tests/test_plugins.py
index 6e31649380..009ef10cc9 100644
--- a/tests/test_plugins.py
+++ b/tests/test_plugins.py
@@ -44,7 +44,7 @@ def test_mode_solver():
         geometry=td.Box(size=(td.inf, 0.5, 0.5)), medium=td.Medium(permittivity=4.0)
     )
     simulation = td.Simulation(size=(2, 2, 2), grid_size=(0.1, 0.1, 0.1), structures=[waveguide])
-    plane = td.Box(center=(0, 0, 0), size=(0, 2, 2))
+    plane = td.Box(center=(0, 0, 0), size=(0, 1, 1))
     ms = ModeSolver(simulation=simulation, plane=plane, freq=td.constants.C_0 / 1.5)
     modes = ms.solve(mode=td.Mode(mode_index=1))
 
diff --git a/tidy3d/components/data.py b/tidy3d/components/data.py
index 8538daa4c4..55799e36f5 100644
--- a/tidy3d/components/data.py
+++ b/tidy3d/components/data.py
@@ -344,7 +344,9 @@ def colocate(self, x, y, z) -> xr.Dataset:
         ----------
         x : np.array
             x coordinates of locations.
+        y : np.array
             y coordinates of locations.
+        z : np.array
             z coordinates of locations.
 
         Returns
diff --git a/tidy3d/components/geometry.py b/tidy3d/components/geometry.py
index 958caaf22e..16749fbadd 100644
--- a/tidy3d/components/geometry.py
+++ b/tidy3d/components/geometry.py
@@ -526,8 +526,8 @@ def from_bounds(cls, rmin: Coordinate, rmax: Coordinate):
         -------
         >>> b = Box.from_bounds(rmin=(-1, -2, -3), rmax=(3, 2, 1))
         """
-        center = tuple((pt_min + pt_max / 2.0) for pt_min, pt_max in zip(rmin, rmax))
-        size = tuple((pt_max - pt_max) for pt_min, pt_max in zip(rmin, rmax))
+        center = tuple((pt_min + pt_max) / 2.0 for pt_min, pt_max in zip(rmin, rmax))
+        size = tuple((pt_max - pt_min) for pt_min, pt_max in zip(rmin, rmax))
         return cls(center=center, size=size)
 
     def intersections(self, x: float = None, y: float = None, z: float = None):
diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py
index 7db3b2a0fd..c9d5687f30 100644
--- a/tidy3d/components/simulation.py
+++ b/tidy3d/components/simulation.py
@@ -142,7 +142,7 @@ class Simulation(Box):  # pylint:disable=too-many-public-methods
     """
     # pylint:enable=line-too-long
 
-    version: str = '0.2.0' # will make this more automated later
+    version: str = "0.2.0"  # will make this more automated later
     grid_size: Tuple[GridSize, GridSize, GridSize]
     medium: MediumType = Medium()
     run_time: pydantic.NonNegativeFloat = 0.0
@@ -825,11 +825,11 @@ def _make_bound_coords_nonuniform(dl, center, size, num_layers):
         bound_coords = np.array([np.sum(dl[:i]) for i in range(len(dl) + 1)])
 
         # shift coords to center at center of simulation along dimension
-        bound_coords = bound_coords - np.sum(dl)/2 + center
+        bound_coords = bound_coords - np.sum(dl) / 2 + center
 
         # chop off any coords outside of simulation bounds
-        bound_min = center - size/2
-        bound_max = center + size/2
+        bound_min = center - size / 2
+        bound_max = center + size / 2
         bound_coords = bound_coords[bound_coords <= bound_max]
         bound_coords = bound_coords[bound_coords >= bound_min]
 
diff --git a/tidy3d/components/source.py b/tidy3d/components/source.py
index d71e3711be..e96c23787b 100644
--- a/tidy3d/components/source.py
+++ b/tidy3d/components/source.py
@@ -12,7 +12,6 @@
 from .geometry import Box
 from .mode import Mode
 from .viz import add_ax_if_none, SourceParams
-from ..log import ValidationError
 from ..constants import inf  # pylint:disable=unused-import
 
 # TODO: change directional source to something signifying its intent is to create a specific field.
@@ -316,7 +315,15 @@ class VolumeSource(Source):
     type: Literal["VolumeSource"] = "VolumeSource"
 
 
-class ModeSource(Source):
+class FieldSource(Source, ABC):
+    """A planar Source defined by the desired E and H fields at a plane. The sources are created
+    such that the propagation is uni-directional."""
+
+    direction: Direction
+    _plane_validator = assert_plane()
+
+
+class ModeSource(FieldSource):
     """Modal profile on finite extent plane
 
     Parameters
@@ -330,9 +337,8 @@ class ModeSource(Source):
     source_time : :class:`GaussianPulse` or :class:`ContinuousWave`
         Specification of time-dependence of source.
     direction : str
-        Specifies the sign of propagation.
-        Must be in ``{'+', '-'}``.
-        Note: propagation occurs along dimension normal to plane.
+        Specifies propagation in the positive or negative direction of the normal axis. Must be in
+        ``{'+', '-'}``.
     mode : :class:`Mode`
         Specification of the mode being injected by source.
     name : str = None
@@ -345,40 +351,11 @@ class ModeSource(Source):
     >>> mode_source = ModeSource(size=(10,10,0), source_time=pulse, mode=mode, direction='-')
     """
 
-    direction: Direction
-    mode: Mode
     type: Literal["ModeSource"] = "ModeSource"
-    _plane_validator = assert_plane()
-
-
-class DirectionalSource(Source, ABC):
-    """A Planar Source with uni-directional propagation."""
-
-    direction: Direction
-    polarization: Polarization
+    mode: Mode
 
-    _plane_validator = assert_plane()
 
-    @pydantic.root_validator(allow_reuse=True)
-    def polarization_is_orthogonal(cls, values):
-        """Ensure the polarization is orthogonal to the propagation direction."""
-        size = values.get("size")
-        polarization = values.get("polarization")
-        assert size is not None
-        assert polarization is not None
-
-        normal_axis_index = size.index(0.0)
-        normal_axis = "xyz"[normal_axis_index]
-        polarization_axis = polarization[-1]
-        if normal_axis == polarization_axis:
-            raise ValidationError(
-                "Directional source must have polarization component orthogonal "
-                "to the normal direction of the plane."
-            )
-        return values
-
-
-class PlaneWave(DirectionalSource):
+class PlaneWave(FieldSource):
     """Uniform distribution on infinite extent plane.
 
     Parameters
@@ -391,14 +368,17 @@ class PlaneWave(DirectionalSource):
         All elements must be non-negative.
     source_time : :class:`GaussianPulse` or :class:`ContinuousWave`
         Specification of time-dependence of source.
-    polarization : str
-        Specifies the direction and type of current component.
-        Must be in ``{'Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz'}``.
-        For example, ``'Ez'`` specifies electric current source polarized along the z-axis.
+    pol_angle : float, optional
+        Specifies the angle between the electric field polarization of the source and the plane
+        defined by the normal axis and the propagation axis (rad). ``pol_angle=0`` (default)
+        specifies P polarization, while ``pol_angle=np.pi/2`` specifies S polarization. At normal
+        incidence when S and P are undefined, ``pol_angle=0`` defines, respectively:
+         - ``Ey`` polarization for propagation along ``x``.
+         - ``Ex`` polarization for propagation along ``y``.
+         - ``Ex`` polarization for propagation along ``z``.
     direction : str
-        Specifies the sign of propagation.
-        Must be in ``{'+', '-'}``.
-        Note: propagation occurs along dimension normal to plane.
+        Specifies propagation in the positive or negative direction of the normal axis. Must be in
+        ``{'+', '-'}``.
     name : str = None
         Optional name for source.
 
@@ -409,9 +389,12 @@ class PlaneWave(DirectionalSource):
     """
 
     type: Literal["PlaneWave"] = "PlaneWave"
+    pol_angle: float = 0
+    # TODO: this is only needed so that the convert path still works. Remove eventually.
+    polarization: Polarization = "Ex"
 
 
-class GaussianBeam(DirectionalSource):
+class GaussianBeam(FieldSource):
     """guassian distribution on finite extent plane
 
     Parameters
@@ -424,26 +407,27 @@ class GaussianBeam(DirectionalSource):
         All elements must be non-negative.
     source_time : :class:`GaussianPulse` or :class:`ContinuousWave`
         Specification of time-dependence of source.
-    polarization : str
-        Specifies the direction and type of current component.
-        Must be in ``{'Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz'}``.
-        For example, ``'Ez'`` specifies electric current source polarized along the z-axis.
     direction : str
-        Specifies the sign of propagation.
-        Must be in ``{'+', '-'}``.
-        Note: propagation occurs along dimension normal to plane.
+        Specifies propagation in the positive or negative direction of the normal axis. Must be in
+        ``{'+', '-'}``.
     waist_radius: float = 1.0
         Radius of the beam at the waist (um).
         Must be positive.
     waist_distance: float = 0.0
         Distance (um) from the beam waist along the propagation direction.
         Must be non-negative.
-    angle_theta: float = 0.0
-        Angle of propagation of the beam with respect to the normal axis (rad).
-    angle_phi: float = 0.0
-        Angle of propagation of the beam with respect to parallel axis (rad).
-    pol_angle: float = 0.0
-        Angle of the polarization with respect to the parallel axis (rad).
+    angle_theta : float, optional
+        Polar angle from the normal axis (rad).
+    angle_phi : float, optional
+        Azimuth angle in the plane orthogonal to the normal axis (rad).
+    pol_angle : float, optional
+        Specifies the angle between the electric field polarization of the source and the plane
+        defined by the normal axis and the propagation axis (rad). ``pol_angle=0`` (default)
+        specifies P polarization, while ``pol_angle=np.pi/2`` specifies S polarization. At normal
+        incidence when S and P are undefined, ``pol_angle=0`` defines, respectively:
+         - ``Ey`` polarization for propagation along ``x``.
+         - ``Ex`` polarization for propagation along ``y``.
+         - ``Ex`` polarization for propagation along ``z``.
     name : str = None
         Optional name for source.
 
@@ -453,7 +437,7 @@ class GaussianBeam(DirectionalSource):
     >>> gauss = GaussianBeam(
     ...     size=(0,3,3),
     ...     source_time=pulse,
-    ...     polarization='Hy',
+    ...     pol_angle=np.pi / 2,
     ...     direction='+',
     ...     waist_radius=1.0)
     """
diff --git a/tidy3d/plugins/mode/dot_product.py b/tidy3d/plugins/mode/dot_product.py
deleted file mode 100644
index 824b9d4313..0000000000
--- a/tidy3d/plugins/mode/dot_product.py
+++ /dev/null
@@ -1,43 +0,0 @@
-# """Dot product between two field distributions in units of power.
-# If fields1 corresponds to a mode normalized such that (fields1, fields1) = 1, then
-# |(fields1, fields2)|^2 is the fraction of the total power carried by fields2 that is
-# specifically carried by the mode given by fields1. """
-
-# import numpy as np
-
-# from ...components import Simulation, Box
-# from .mode_solver import ModeInfo
-
-# def dot_product(sim: Simulation, plane: Box, mode1: ModeInfo, mode2: ModeInfo) -> float:
-#     """Dot product between two modes.
-
-#     Parameters
-#     ----------
-#     sim : Simulation
-#         Simulation in which the modes were computed.
-#     plane : Box
-#         The plane at which the modes were computed.
-#     mode1 : ModeInfo
-#         Data structure of the first mode.
-#     mode2 : ModeInfo
-#         Data structure of the second mode.
-
-#     Returns
-#     -------
-#     float
-#         The overlap integral between the two modes.
-#     """
-
-#     E1 = fields1[0][:2, :, :]
-#     H1 = fields1[1][:2, :, :]
-#     E2 = fields2[0][:2, :, :]
-#     H2 = fields2[1][:2, :, :]
-
-#     dl1 = coords[0][1:] - coords[0][:-1]
-#     dl2 = coords[1][1:] - coords[1][:-1]
-#     dA = np.outer(dl1, dl2)
-#     dV = dA * (coords[2][1] - coords[2][0]) ** (1 / 4)
-
-#     cross = np.cross(np.conj(E1), H2, axis=0) + np.cross(E2, np.conj(H1), axis=0)
-
-#     return 1 / 4 * np.sum(cross * dA)
diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py
index 094c0c167e..c6f3f3be39 100644
--- a/tidy3d/plugins/mode/mode_solver.py
+++ b/tidy3d/plugins/mode/mode_solver.py
@@ -10,7 +10,6 @@
 from ...components import Box
 from ...components import Simulation
 from ...components import Mode
-from ...components import FieldData, ScalarFieldData
 from ...components import ModeMonitor
 from ...components import ModeSource, GaussianPulse
 from ...components.types import Direction
@@ -50,10 +49,9 @@
 @dataclass
 class ModeInfo:
     """stores information about a (solved) mode.
-
     Attributes
     ----------
-    field_Data: FieldData
+    field_data: xr.Dataset
         Contains information about the fields of the modal profile.
     mode: Mode
         Specifications of the mode.
@@ -202,26 +200,25 @@ def rotate_field_coords(field_array):
         for field_name, field in fields.items():
             plane_grid = self.simulation.discretize(self.plane)
             plane_coords = plane_grid[field_name]
-
-            data_dict[field_name] = ScalarFieldData(
-                x=plane_coords.x,
-                y=plane_coords.y,
-                z=plane_coords.z,
-                f=np.array([self.freq]),
-                values=field,
-            )
+            coords = {
+                "x": plane_coords.x,
+                "y": plane_coords.y,
+                "z": plane_coords.z,
+                "f": np.array([self.freq]),
+            }
+            data_dict[field_name] = xr.DataArray(field, coords=coords)
 
         n_eff_complex = n_eff_complex[mode.mode_index]
 
-        field_data = FieldData(data_dict=data_dict)
-
-        return ModeInfo(
-            field_data=field_data.data,
+        mode_info = ModeInfo(
+            field_data=xr.Dataset(data_dict),
             mode=mode,
             n_eff=n_eff_complex.real,
             k_eff=n_eff_complex.imag,
         )
 
+        return mode_info
+
     def make_source(self, mode: Mode, fwidth: float, direction: Direction) -> ModeSource:
         """Creates ``ModeSource`` from a Mode + additional specifications.