diff --git a/.github/workflows/conda-forge.yml b/.github/workflows/conda-forge.yml
deleted file mode 100644
index 121cd6db94d..00000000000
--- a/.github/workflows/conda-forge.yml
+++ /dev/null
@@ -1,73 +0,0 @@
-#  For more information about TARDIS pipelines, please refer to:
-#
-#    https://tardis-sn.github.io/tardis/contributing/development/continuous_integration.html
-
-name: conda-forge
-
-on:
-  schedule:
-    - cron: '0 0 * * 1'
-
-  workflow_dispatch:  # manual trigger
-
-defaults:
-  run:
-    shell: bash -l {0}
-
-jobs:
-  build:
-    strategy:
-      fail-fast: false
-      matrix:
-        pkgman: [conda, mamba]
-        os: [ubuntu-latest, macos-latest]
-
-        include:
-        - os: ubuntu-latest
-          platform: linux-64
-
-        - os: macos-latest
-          platform: osx-64
-
-    runs-on: ${{ matrix.os }}
-    name: ${{ matrix.platform }} (${{ matrix.pkgman }})
-    steps:
-      - name: Download lockfile
-        run: wget -q https://github.com/tardis-sn/tardis/releases/latest/download/conda-${{ matrix.platform }}.lock
-
-      - name: Setup environment
-        uses: conda-incubator/setup-miniconda@v2
-        with:
-            miniforge-variant: Mambaforge
-            miniforge-version: latest
-            environment-file: conda-${{ matrix.platform }}.lock
-            activate-environment: tardis
-            use-mamba: true
-
-      - name: Install conda-forge package
-        run: ${{ matrix.pkgman }} install tardis-sn -c conda-forge
-
-  slack-notify:
-    if: ${{ always() }}
-    needs: build
-    runs-on: ubuntu-latest
-    steps:
-      - name: Send custom JSON data to Slack workflow
-        uses: slackapi/slack-github-action@v1
-        with:
-          payload: |
-            {
-              "text": "conda-forge package installation: ${{ needs.build.result }}\n${{ github.event.pull_request.html_url || github.event.head_commit.url }}",
-              "blocks": [
-                {
-                  "type": "section",
-                  "text": {
-                    "type": "mrkdwn",
-                    "text": "conda-forge package installation: ${{ needs.build.result }}\n${{ github.event.pull_request.html_url || github.event.head_commit.url }}"
-                  }
-                }
-              ]
-            }
-        env:
-          SLACK_WEBHOOK_URL: ${{ secrets.SLACK_WEBHOOK_URL }}
-          SLACK_WEBHOOK_TYPE: INCOMING_WEBHOOK
diff --git a/.zenodo.json b/.zenodo.json
index 06827c53e3b..8720d4b208a 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -116,14 +116,17 @@
       "name": "Sondhi, Dhruv"
     },
     {
-      "name": "Barbosa, Talytha"
+      "name": "O'Brien, Jack"
     },
     {
-      "name": "O'Brien, Jack"
+      "name": "Barbosa, Talytha"
     },
     {
       "name": "Yu, Jenny"
     },
+    {
+      "name": "Shields, Joshua"
+    },
     {
       "name": "Patel, Maryam"
     },
@@ -142,20 +145,20 @@
       "affiliation": "University of Oxford"
     },
     {
-      "name": "Shields, Joshua"
+      "name": "Singh, Shreyas"
     },
     {
       "name": "Savel, Arjun"
     },
-    {
-      "name": "Singh, Shreyas"
-    },
     {
       "name": "Reinecke, Martin"
     },
     {
       "name": "Eweis, Youssef"
     },
+    {
+      "name": "Holas, Alexander"
+    },
     {
       "name": "Bylund, Tomas"
     },
@@ -165,21 +168,18 @@
     {
       "name": "Bentil, Laud"
     },
-    {
-      "name": "Holas, Alexander"
-    },
     {
       "name": "Eguren, Jordi",
       "orcid": "0000-0002-2328-8030"
     },
     {
-      "name": "Alam, Arib"
+      "name": "Bartnik, Matthew"
     },
     {
-      "name": "Kumar, Ansh"
+      "name": "Alam, Arib"
     },
     {
-      "name": "Bartnik, Matthew"
+      "name": "Kumar, Ansh"
     },
     {
       "name": "Varma Buddaraju, Rohith"
@@ -187,15 +187,11 @@
     {
       "name": "Magee, Mark"
     },
-    {
-      "name": "Livneh, Ran"
-    },
     {
       "name": "Kambham, Satwik"
     },
     {
-      "name": "Mishra, Sashank",
-      "orcid": "0000-0001-8302-1584"
+      "name": "Livneh, Ran"
     },
     {
       "name": "Rajagopalan, Srinath"
@@ -204,37 +200,38 @@
       "name": "Daksh, Ayushi"
     },
     {
-      "name": "Bhakar, Jayant"
+      "name": "Mishra, Sashank",
+      "orcid": "0000-0001-8302-1584"
     },
     {
-      "name": "Reichenbach, John"
+      "name": "Dutta, Anirban"
     },
     {
       "name": "Floers, Andreas"
     },
     {
-      "name": "Actions, GitHub"
+      "name": "Jain, Rinkle"
     },
     {
-      "name": "Jain, Rinkle"
+      "name": "Reichenbach, John"
     },
     {
-      "name": "Singh, Sourav"
+      "name": "Actions, GitHub"
     },
     {
-      "name": "Brar, Antreev"
+      "name": "Bhakar, Jayant"
     },
     {
-      "name": "Dutta, Anirban"
+      "name": "Brar, Antreev"
     },
     {
       "name": "Chaumal, Aarya"
     },
     {
-      "name": "Matsumura, Yuki"
+      "name": "Singh, Sourav"
     },
     {
-      "name": "Patidar, Abhishek"
+      "name": "Selsing, Jonatan"
     },
     {
       "name": "Kowalski, Nathan"
@@ -246,16 +243,19 @@
       "name": "Gupta, Harshul"
     },
     {
-      "name": "Selsing, Jonatan"
+      "name": "Talegaonkar, Chinmay"
     },
     {
-      "name": "Talegaonkar, Chinmay"
+      "name": "Matsumura, Yuki"
     },
     {
       "name": "Sofiatti, Caroline"
     },
     {
-      "name": "Prasad, Shilpi"
+      "name": "Patidar, Abhishek"
+    },
+    {
+      "name": "Venkat, Shashank"
     },
     {
       "name": "Buchner, Johannes"
@@ -263,6 +263,9 @@
     {
       "name": "Yap, Kevin"
     },
+    {
+      "name": "Martinez, Laureano"
+    },
     {
       "name": "Truong, Le"
     },
@@ -276,22 +279,22 @@
       "name": "Sarafina, Nance"
     },
     {
-      "name": "Patra, Nilesh"
+      "name": "Volodin, Dmitry"
     },
     {
-      "name": "Volodin, Dmitry"
+      "name": "Patra, Nilesh"
     },
     {
-      "name": "Singh Rathore, Parikshit"
+      "name": "Patel, Pratik"
     },
     {
       "name": "Sharma, Sampark"
     },
     {
-      "name": "Venkat, Shashank"
+      "name": "Lu, Jing"
     },
     {
-      "name": "Lu, Jing"
+      "name": "Prasad, Shilpi"
     },
     {
       "name": "Gupta, Suyash"
@@ -306,7 +309,7 @@
       "name": "Aggarwal, Yash"
     },
     {
-      "name": "Patel, Pratik"
+      "name": "Singh Rathore, Parikshit"
     },
     {
       "name": "Dasgupta, Debajyoti"
@@ -315,19 +318,16 @@
       "name": "PATIDAR, ABHISHEK"
     },
     {
-      "name": "Martinez, Laureano"
+      "name": "Nayak U, Ashwin"
     },
     {
       "name": "Kolliboyina, Chaitanya"
     },
     {
-      "name": "Nayak U, Ashwin"
+      "name": "Kumar, Atul"
     },
     {
       "name": "Kharkar, Atharwa"
-    },
-    {
-      "name": "Kumar, Atul"
     }
   ]
 }
\ No newline at end of file
diff --git a/docs/io/optional/callback_example.ipynb b/docs/io/optional/callback_example.ipynb
index 0b13d5bb8cd..04fa294f37e 100644
--- a/docs/io/optional/callback_example.ipynb
+++ b/docs/io/optional/callback_example.ipynb
@@ -78,7 +78,7 @@
    "source": [
     "def append_num_emitted_to_list(sim, lst):\n",
     "    if sim.iterations_executed < sim.iterations:\n",
-    "        num_emitted_packets = len(sim.transport.emitted_packet_nu)\n",
+    "        num_emitted_packets = len(sim.transport.transport_state.emitted_packet_nu)\n",
     "        lst.append(num_emitted_packets)"
    ]
   },
diff --git a/docs/io/optional/custom_source.ipynb b/docs/io/optional/custom_source.ipynb
index 101df3df447..767421a792d 100644
--- a/docs/io/optional/custom_source.ipynb
+++ b/docs/io/optional/custom_source.ipynb
@@ -42,6 +42,9 @@
     "from tardis import constants as const\n",
     "from astropy import units as u\n",
     "from tardis.montecarlo.packet_source import BlackBodySimpleSource\n",
+    "from tardis.montecarlo.montecarlo_numba.packet_collections import (\n",
+    "    PacketCollection,\n",
+    ")\n",
     "from tardis import run_tardis\n",
     "import matplotlib.pyplot as plt\n",
     "from tardis.io.atom_data import download_atom_data"
@@ -87,7 +90,7 @@
     "        self.truncation_wavelength = truncation_wavelength\n",
     "        super().__init__(**kwargs)\n",
     "\n",
-    "    def create_packets(self, no_of_packets, drawing_sample_size=None):\n",
+    "    def create_packets(self, no_of_packets, drawing_sample_size=None, seed_offset=0, *args, **kwargs):\n",
     "        \"\"\"\n",
     "        Packet source that generates a truncated Blackbody source.\n",
     "\n",
@@ -108,12 +111,17 @@
     "            Packet energies\n",
     "        \"\"\"\n",
     "\n",
+    "        self._reseed(self.base_seed + seed_offset)\n",
+    "        packet_seeds = self.rng.choice(\n",
+    "            self.MAX_SEED_VAL, no_of_packets, replace=True\n",
+    "        )\n",
+    "\n",
     "        # Makes uniform array of packet radii from blackbody source\n",
-    "        radii = self.create_packet_radii(no_of_packets)\n",
+    "        radii = self.create_packet_radii(no_of_packets, *args, **kwargs)\n",
     "\n",
     "        # Use mus and energies from normal blackbody source.\n",
-    "        mus = self.create_packet_mus(no_of_packets)\n",
-    "        energies = self.create_packet_energies(no_of_packets)\n",
+    "        mus = self.create_packet_mus(no_of_packets, *args, **kwargs)\n",
+    "        energies = self.create_packet_energies(no_of_packets, *args, **kwargs)\n",
     "\n",
     "        # If not specified, draw 2 times as many packets and reject any beyond no_of_packets.\n",
     "        if drawing_sample_size is None:\n",
@@ -128,18 +136,22 @@
     "\n",
     "        # Draw nus from blackbody distribution and reject based on truncation_frequency.\n",
     "        # If more nus.shape[0] > no_of_packets use only the first no_of_packets.\n",
-    "        nus = self.create_packet_nus(drawing_sample_size)\n",
+    "        nus = self.create_packet_nus(drawing_sample_size, *args, **kwargs)\n",
     "        nus = nus[nus < truncation_frequency][:no_of_packets]\n",
     "\n",
     "        # Only required if the truncation wavelength is too big compared to the maximum\n",
     "        # of the blackbody distribution. Keep sampling until nus.shape[0] > no_of_packets.\n",
     "        while nus.shape[0] < no_of_packets:\n",
-    "            additional_nus = self.create_packet_nus(drawing_sample_size)\n",
+    "            additional_nus = self.create_packet_nus(drawing_sample_size, *args, **kwargs)\n",
     "            mask = additional_nus < truncation_frequency\n",
     "            additional_nus = additional_nus[mask][:no_of_packets]\n",
     "            nus = np.hstack([nus, additional_nus])[:no_of_packets]\n",
     "\n",
-    "        return radii, nus, mus, energies"
+    "        radiation_field_luminosity = (\n",
+    "            self.calculate_radfield_luminosity().to(u.erg / u.s).value\n",
+    "        )\n",
+    "\n",
+    "        return PacketCollection(radii, nus, mus, energies, packet_seeds, radiation_field_luminosity)"
    ]
   },
   {
@@ -179,11 +191,11 @@
    "outputs": [],
    "source": [
     "%matplotlib inline\n",
-    "plt.plot(mdl.transport.spectrum_virtual.wavelength,\n",
-    "         mdl.transport.spectrum_virtual.luminosity_density_lambda,\n",
+    "plt.plot(mdl.transport.transport_state.spectrum_virtual.wavelength,\n",
+    "         mdl.transport.transport_state.spectrum_virtual.luminosity_density_lambda,\n",
     "         color='red', label='truncated blackbody (custom packet source)')\n",
-    "plt.plot(mdl_norm.transport.spectrum_virtual.wavelength,\n",
-    "         mdl_norm.transport.spectrum_virtual.luminosity_density_lambda,\n",
+    "plt.plot(mdl_norm.transport.transport_state.spectrum_virtual.wavelength,\n",
+    "         mdl_norm.transport.transport_state.spectrum_virtual.luminosity_density_lambda,\n",
     "         color='blue', label='normal blackbody (default packet source)')\n",
     "plt.xlabel('$\\lambda [\\AA]$')\n",
     "plt.ylabel('$L_\\lambda$ [erg/s/$\\AA$]')\n",
@@ -208,7 +220,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.10"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,
diff --git a/docs/io/output/rpacket_tracking.ipynb b/docs/io/output/rpacket_tracking.ipynb
index fc4b1dc002e..b2a1127f00f 100644
--- a/docs/io/output/rpacket_tracking.ipynb
+++ b/docs/io/output/rpacket_tracking.ipynb
@@ -202,7 +202,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "type(sim.transport.rpacket_tracker)"
+    "type(sim.transport.transport_state.rpacket_tracker)"
    ]
   },
   {
@@ -243,7 +243,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "len(sim.transport.rpacket_tracker)"
+    "len(sim.transport.transport_state.rpacket_tracker)"
    ]
   },
   {
@@ -281,7 +281,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker[10].index"
+    "sim.transport.transport_state.rpacket_tracker[10].index"
    ]
   },
   {
@@ -300,7 +300,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker[10].seed"
+    "sim.transport.transport_state.rpacket_tracker[10].seed"
    ]
   },
   {
@@ -319,7 +319,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker[10].status"
+    "sim.transport.transport_state.rpacket_tracker[10].status"
    ]
   },
   {
@@ -347,7 +347,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "len(sim.transport.rpacket_tracker[10].shell_id)"
+    "len(sim.transport.transport_state.rpacket_tracker[10].shell_id)"
    ]
   },
   {
@@ -403,7 +403,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "type(sim.transport.rpacket_tracker_df)"
+    "type(sim.transport.transport_state.rpacket_tracker_df)"
    ]
   },
   {
@@ -461,7 +461,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker_df"
+    "sim.transport.transport_state.rpacket_tracker_df"
    ]
   },
   {
@@ -490,7 +490,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker_df.loc[10]"
+    "sim.transport.transport_state.rpacket_tracker_df.loc[10]"
    ]
   },
   {
@@ -509,7 +509,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker_df.loc[10][\"energy\"]"
+    "sim.transport.transport_state.rpacket_tracker_df.loc[10][\"energy\"]"
    ]
   },
   {
@@ -528,7 +528,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker_df.loc[10,5][\"energy\"]"
+    "sim.transport.transport_state.rpacket_tracker_df.loc[10,5][\"energy\"]"
    ]
   },
   {
@@ -538,7 +538,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "sim.transport.rpacket_tracker_df.loc[10][\"energy\"][5]"
+    "sim.transport.transport_state.rpacket_tracker_df.loc[10][\"energy\"][5]"
    ]
   },
   {
@@ -567,7 +567,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.13"
+   "version": "3.11.5"
   },
   "vscode": {
    "interpreter": {
diff --git a/docs/io/visualization/montecarlo_packet_visualization.ipynb b/docs/io/visualization/montecarlo_packet_visualization.ipynb
index fe46b67b745..1c15aac1657 100644
--- a/docs/io/visualization/montecarlo_packet_visualization.ipynb
+++ b/docs/io/visualization/montecarlo_packet_visualization.ipynb
@@ -89,7 +89,7 @@
    "outputs": [],
    "source": [
     "# accessing the rpacket_tracker dataframe\n",
-    "sim.transport.rpacket_tracker_df"
+    "sim.transport.transport_state.rpacket_tracker_df"
    ]
   },
   {
@@ -191,11 +191,11 @@
    "source": [
     "# function for getting coordinates of all packets\n",
     "def get_coordinates_multiple_packets_ints(r_packet_tracker,time):\n",
-    "    thetas = np.linspace(0, 2*math.pi,sim.transport.rpacket_tracker_df.index[-1][0]+2)\n",
+    "    thetas = np.linspace(0, 2*math.pi,sim.transport.transport_state.rpacket_tracker_df.index[-1][0]+2)\n",
     "    x = []\n",
     "    y = []\n",
     "    inters = []\n",
-    "    for i in range(sim.transport.rpacket_tracker_df.index[-1][0]+1):\n",
+    "    for i in range(sim.transport.transport_state.rpacket_tracker_df.index[-1][0]+1):\n",
     "        xs,ys,ints = get_x_y_ints_with_theta_init(r_packet_tracker.loc[i][\"r\"],r_packet_tracker.loc[i][\"mu\"],time,r_packet_tracker.loc[i][\"interaction_type\"],'n',thetas[i])\n",
     "        x.append(xs)\n",
     "        y.append(ys)\n",
@@ -236,7 +236,7 @@
     "# animated plot\n",
     "\n",
     "v_shells = sim.simulation_state.radius.value * 1e-5 / sim.simulation_state.time_explosion.value\n",
-    "xs,ys,ints = get_coordinates_multiple_packets_ints(sim.transport.rpacket_tracker_df,sim.simulation_state.time_explosion.value)\n",
+    "xs,ys,ints = get_coordinates_multiple_packets_ints(sim.transport.transport_state.rpacket_tracker_df,sim.simulation_state.time_explosion.value)\n",
     "xs,ys,ints,max_size = get_equal_array_size(xs,ys,ints)\n",
     "\n",
     "fig = go.Figure()\n",
@@ -273,7 +273,7 @@
     "\n",
     "#Add the packet trajectory\n",
     "\n",
-    "df = sim.transport.rpacket_tracker_df\n",
+    "df = sim.transport.transport_state.rpacket_tracker_df\n",
     "interaction_from_num = {0: \"No Interaction\", 1: \"EScattering\", 2: \"Line\"}\n",
     "interaction_color_from_num = {0: \"darkslategrey\", 1: \"#3366FF\", 2: \"#FF3300\"}\n",
     "interaction_opacity_from_num = {0: 0, 1: 1, 2: 1}\n",
diff --git a/docs/physics/montecarlo/initialization.ipynb b/docs/physics/montecarlo/initialization.ipynb
index f53662dee62..2ab262e6da1 100644
--- a/docs/physics/montecarlo/initialization.ipynb
+++ b/docs/physics/montecarlo/initialization.ipynb
@@ -79,6 +79,9 @@
    "source": [
     "import numpy as np\n",
     "from tardis.montecarlo.packet_source import BlackBodySimpleSource\n",
+    "from tardis.montecarlo.montecarlo_numba.packet_collections import (\n",
+    "    PacketCollection,\n",
+    ")\n",
     "from astropy import units as u\n",
     "from tardis import constants as const\n",
     "import matplotlib.pyplot as plt"
@@ -164,24 +167,21 @@
     "# We define our packet source\n",
     "packet_source = BlackBodySimpleSource(base_seed=base_seed)\n",
     "\n",
-    "# Create separate packet seeds to sample each packet's montecarlo trajectory\n",
-    "packet_seeds = packet_source.create_packet_seeds(n_packets, seed_offset)\n",
-    "\n",
     "# Set radii and temperature from model\n",
     "packet_source.radius = r_boundary_inner\n",
-    "packet_source.temperature = temperature_inner.value\n",
+    "packet_source.temperature = temperature_inner\n",
     "\n",
     "# Create packets\n",
-    "radii, nus, mus, energies = packet_source.create_packets(n_packets)\n",
+    "packet_collection = packet_source.create_packets(n_packets)\n",
     "\n",
     "# Sets the energies in units of ergs\n",
-    "energies *= u.erg\n",
+    "packet_collection.initial_energies *= u.erg\n",
     "\n",
     "# Sets the frequencies in units of Hz\n",
-    "nus *= u.Hz\n",
+    "packet_collection.initial_nus *= u.Hz\n",
     "\n",
-    "print(\"Energies:\", energies)\n",
-    "print(\"Radii:\", radii)"
+    "print(\"Energies:\", packet_collection.initial_energies)\n",
+    "print(\"Radii:\", packet_collection.initial_radii)"
    ]
   },
   {
@@ -205,7 +205,7 @@
     "print(\"Time of simulation:\", t_simulation)\n",
     "\n",
     "# Array of luminosity contribution by each packet\n",
-    "lumin_per_packet = energies / t_simulation\n",
+    "lumin_per_packet = packet_collection.initial_energies / t_simulation\n",
     "print(\"Luminosity per packet:\", lumin_per_packet)"
    ]
   },
@@ -263,15 +263,15 @@
    "source": [
     "# We set important quantites for making our histogram\n",
     "bins = 200\n",
-    "nus_planck = np.linspace(min(nus), max(nus), bins)\n",
+    "nus_planck = np.linspace(min(packet_collection.initial_nus), max(packet_collection.initial_nus), bins).value\n",
     "bin_width = nus_planck[1] - nus_planck[0]\n",
     "\n",
     "# In the histogram plot below, the weights argument is used\n",
     "# to make sure our plotted spectrum has the correct y-axis scale\n",
-    "plt.hist(nus.value, bins=bins, weights=lumin_per_packet / bin_width)\n",
+    "plt.hist(packet_collection.initial_nus.value, bins=bins, weights=lumin_per_packet / bin_width)\n",
     "\n",
     "# We plot the planck function for comparison\n",
-    "plt.plot(nus_planck, planck_function(nus_planck))\n",
+    "plt.plot(nus_planck * u.Hz, planck_function(nus_planck * u.Hz))\n",
     "\n",
     "plt.xlabel(\"Frequency (Hz)\")\n",
     "plt.ylabel(\"Luminosity density w.r.t. frequency (erg/s/Hz)\")\n",
@@ -296,7 +296,7 @@
    "source": [
     "x = np.linspace(0, 1, 1000)\n",
     "\n",
-    "plt.hist(mus, bins=bins, density=True)\n",
+    "plt.hist(packet_collection.initial_mus, bins=bins, density=True)\n",
     "plt.plot(x, 2 * x)\n",
     "plt.xlabel(\"Propagation direction\")\n",
     "plt.ylabel(\"Probability density\")\n",
@@ -312,7 +312,7 @@
    "source": [
     "thetas = np.linspace(0, np.pi / 2, 1000)\n",
     "\n",
-    "plt.hist(np.arccos(mus), bins=bins, density=True)\n",
+    "plt.hist(np.arccos(packet_collection.initial_mus), bins=bins, density=True)\n",
     "plt.plot(thetas, np.sin(2 * thetas))\n",
     "plt.xlabel(\"Angle with normal (rad)\")\n",
     "plt.ylabel(\"Probability density\")\n",
@@ -348,7 +348,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.13"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,
diff --git a/docs/physics/spectrum/basic.ipynb b/docs/physics/spectrum/basic.ipynb
index d0d61571ea7..2ac6c24799f 100644
--- a/docs/physics/spectrum/basic.ipynb
+++ b/docs/physics/spectrum/basic.ipynb
@@ -47,7 +47,7 @@
    "source": [
     "from tardis.io.configuration.config_reader import Configuration\n",
     "from tardis.simulation import Simulation\n",
-    "from tardis.montecarlo import TARDISSpectrum\n",
+    "from tardis.montecarlo.spectrum import TARDISSpectrum\n",
     "from tardis.io.atom_data.util import download_atom_data\n",
     "from astropy import units as u\n",
     "import numpy as np\n",
@@ -118,7 +118,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "nus = sim.transport.output_nu\n",
+    "nus = sim.transport.transport_state.output_nu\n",
     "nus"
    ]
   },
@@ -129,7 +129,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "energies = sim.transport.output_energy\n",
+    "energies = sim.transport.transport_state.output_energy\n",
     "energies"
    ]
   },
@@ -183,7 +183,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "luminosities = energies / sim.transport.time_of_simulation\n",
+    "luminosities = energies / sim.transport.transport_state.time_of_simulation\n",
     "luminosities"
    ]
   },
@@ -203,7 +203,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "emitted_mask = sim.transport.emitted_packet_mask\n",
+    "emitted_mask = sim.transport.transport_state.emitted_packet_mask\n",
     "emitted_mask"
    ]
   },
@@ -460,7 +460,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.12"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,
diff --git a/docs/physics/update_and_conv/update_and_conv.ipynb b/docs/physics/update_and_conv/update_and_conv.ipynb
index 971249e6a75..b0c8da0d24c 100644
--- a/docs/physics/update_and_conv/update_and_conv.ipynb
+++ b/docs/physics/update_and_conv/update_and_conv.ipynb
@@ -311,7 +311,7 @@
    },
    "outputs": [],
    "source": [
-    "j_estimator = transport.j_estimator * (u.erg * u.cm) \n",
+    "j_estimator = transport.transport_state.estimators.j_estimator * (u.erg * u.cm) \n",
     "j_estimator"
    ]
   },
@@ -324,7 +324,7 @@
    },
    "outputs": [],
    "source": [
-    "nu_bar_estimator = transport.nu_bar_estimator * (u.erg * u.cm * u.Hz)\n",
+    "nu_bar_estimator = transport.transport_state.estimators.nu_bar_estimator * (u.erg * u.cm * u.Hz)\n",
     "nu_bar_estimator"
    ]
   },
@@ -345,7 +345,7 @@
    "outputs": [],
    "source": [
     "V = simulation_state.volume\n",
-    "Delta_t = transport.calculate_time_of_simulation(simulation_state)\n",
+    "Delta_t = transport.transport_state.time_of_simulation\n",
     "prefactor = 1 / (4 * np.pi * V * Delta_t)\n",
     "J = prefactor * j_estimator\n",
     "J"
@@ -462,7 +462,7 @@
     "#nu_lower = tardis_config.supernova.luminosity_wavelength_end.to(u.Hz, u.spectral)\n",
     "#nu_upper = tardis_config.supernova.luminosity_wavelength_start.to(u.Hz, u.spectral)\n",
     "\n",
-    "L_output = transport.calculate_emitted_luminosity(0,np.inf)\n",
+    "L_output = transport.transport_state.calculate_emitted_luminosity(0,np.inf)\n",
     "L_output"
    ]
   },
@@ -570,6 +570,14 @@
    "source": [
     "plasma.electron_densities"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bd488259",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -588,7 +596,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.13"
+   "version": "3.11.5"
   },
   "vscode": {
    "interpreter": {
diff --git a/docs/quickstart.ipynb b/docs/quickstart.ipynb
index b694d9dcca1..33e69d0d773 100644
--- a/docs/quickstart.ipynb
+++ b/docs/quickstart.ipynb
@@ -145,9 +145,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "spectrum = sim.transport.spectrum\n",
-    "spectrum_virtual = sim.transport.spectrum_virtual\n",
-    "spectrum_integrated = sim.transport.spectrum_integrated"
+    "spectrum = sim.transport.transport_state.spectrum\n",
+    "spectrum_virtual = sim.transport.transport_state.spectrum_virtual\n",
+    "spectrum_integrated = sim.transport.transport_state.spectrum_integrated"
    ]
   },
   {
@@ -191,7 +191,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.17"
+   "version": "3.11.5"
   }
  },
  "nbformat": 4,
diff --git a/tardis/conftest.py b/tardis/conftest.py
index 5c45a307db9..42530086111 100644
--- a/tardis/conftest.py
+++ b/tardis/conftest.py
@@ -232,3 +232,14 @@ def simulation_verysimple(config_verysimple, atomic_dataset):
     sim = Simulation.from_config(config_verysimple, atom_data=atomic_data)
     sim.iterate(4000)
     return sim
+
+
+@pytest.fixture(scope="session")
+def simulation_verysimple_vpacket_tracking(config_verysimple, atomic_dataset):
+    atomic_data = deepcopy(atomic_dataset)
+    sim = Simulation.from_config(
+        config_verysimple, atom_data=atomic_data, virtual_packet_logging=True
+    )
+    sim.last_no_of_packets = 4000
+    sim.run_final()
+    return sim
diff --git a/tardis/io/model_reader.py b/tardis/io/model_reader.py
index 5390cd24e7d..67beb541ab2 100644
--- a/tardis/io/model_reader.py
+++ b/tardis/io/model_reader.py
@@ -12,7 +12,7 @@
 from radioactivedecay.utils import Z_DICT, elem_to_Z
 
 from tardis.io.configuration.config_reader import ConfigurationNameSpace
-from tardis.montecarlo.base import MontecarloTransportSolver
+from tardis.montecarlo.base import MonteCarloTransportSolver
 from tardis.montecarlo.packet_source import (
     BlackBodySimpleSource,
     BlackBodySimpleSourceRelativistic,
@@ -726,7 +726,7 @@ def transport_from_hdf(fname):
     )
 
     # Creating a transport object and storing data
-    new_transport = MontecarloTransportSolver(
+    new_transport = MonteCarloTransportSolver(
         spectrum_frequency=d["spectrum_frequency_cgs"],
         virtual_spectrum_spawn_range=d["virtual_spectrum_spawn_range"],
         disable_electron_scattering=d["disable_electron_scattering"],
@@ -739,7 +739,7 @@ def transport_from_hdf(fname):
         spectrum_method=d["spectrum_method"],
         packet_source=d["packet_source"],
         nthreads=d["nthreads"],
-        enable_vpacket_tracking=d["virt_logging"],
+        enable_virtual_packet_logging=d["virt_logging"],
         use_gpu=d["use_gpu"],
     )
 
diff --git a/tardis/model/parse_input.py b/tardis/model/parse_input.py
index 63b3658ec30..ceddfa04984 100644
--- a/tardis/model/parse_input.py
+++ b/tardis/model/parse_input.py
@@ -20,7 +20,10 @@
 from tardis.model.radiation_field_state import (
     DiluteBlackBodyRadiationFieldState,
 )
-from tardis.montecarlo.packet_source import BlackBodySimpleSource
+from tardis.montecarlo.packet_source import (
+    BlackBodySimpleSource,
+    BlackBodySimpleSourceRelativistic,
+)
 from tardis.util.base import quantity_linspace
 
 logger = logging.getLogger(__name__)
@@ -574,9 +577,9 @@ def parse_radiation_field_state(
     return DiluteBlackBodyRadiationFieldState(t_radiative, dilution_factor)
 
 
-def parse_packet_source(config, geometry):
+def initialize_packet_source(config, geometry, packet_source):
     """
-    Parse the packet source based on the given configuration and geometry.
+    Initialize the packet source based on config and geometry
 
     Parameters
     ----------
@@ -584,26 +587,19 @@ def parse_packet_source(config, geometry):
         The configuration object containing the supernova and plasma settings.
     geometry : Geometry
         The geometry object containing the inner radius information.
+    packet_source : BasePacketSource
+        The packet source object based on the configuration and geometry.
 
     Returns
     -------
-    packet_source : BlackBodySimpleSource
+    packet_source : BasePacketSource
         The packet source object based on the configuration and geometry.
 
     Raises
     ------
     ValueError
         If both t_inner and luminosity_requested are None.
-
     """
-    if config.montecarlo.enable_full_relativity:
-        packet_source = BlackBodySimpleSourceRelativistic(
-            base_seed=config.montecarlo.seed,
-            time_explosion=config.supernova.time_explosion,
-        )
-    else:
-        packet_source = BlackBodySimpleSource(base_seed=config.montecarlo.seed)
-
     luminosity_requested = config.supernova.luminosity_requested
     if config.plasma.initial_t_inner > 0.0 * u.K:
         packet_source.radius = geometry.r_inner[0]
@@ -620,6 +616,32 @@ def parse_packet_source(config, geometry):
         )
     return packet_source
 
+def parse_packet_source(config, geometry):
+    """
+    Parse the packet source based on the given configuration and geometry.
+
+    Parameters
+    ----------
+    config : Config
+        The configuration object containing the supernova and plasma settings.
+    geometry : Geometry
+        The geometry object containing the inner radius information.
+
+    Returns
+    -------
+    packet_source : BlackBodySimpleSource
+        The packet source object based on the configuration and geometry.
+    """
+    if config.montecarlo.enable_full_relativity:
+        packet_source = BlackBodySimpleSourceRelativistic(
+            base_seed=config.montecarlo.seed,
+            time_explosion=config.supernova.time_explosion,
+        )
+    else:
+        packet_source = BlackBodySimpleSource(base_seed=config.montecarlo.seed)
+
+    return initialize_packet_source(config, geometry, packet_source)
+
 
 def parse_csvy_radiation_field_state(
     config, csvy_model_config, csvy_model_data, geometry, packet_source
diff --git a/tardis/montecarlo/base.py b/tardis/montecarlo/base.py
index 0956e2c3245..6820c61b98d 100644
--- a/tardis/montecarlo/base.py
+++ b/tardis/montecarlo/base.py
@@ -39,26 +39,14 @@
 
 
 # TODO: refactor this into more parts
-class MontecarloTransportSolver(HDFWriterMixin):
+class MonteCarloTransportSolver(HDFWriterMixin):
     """
-    This class is designed as an interface between the Python part and the
-    montecarlo C-part
+    This class modifies the MonteCarloTransportState to solve the radiative
+    transfer problem.
     """
 
     hdf_properties = ["transport_state"]
 
-    vpacket_hdf_properties = [
-        "virt_packet_nus",
-        "virt_packet_energies",
-        "virt_packet_initial_rs",
-        "virt_packet_initial_mus",
-        "virt_packet_last_interaction_in_nu",
-        "virt_packet_last_interaction_type",
-        "virt_packet_last_line_interaction_in_id",
-        "virt_packet_last_line_interaction_out_id",
-        "virt_packet_last_line_interaction_shell_id",
-    ]
-
     hdf_name = "transport"
 
     def __init__(
@@ -74,7 +62,7 @@ def __init__(
         v_packet_settings,
         spectrum_method,
         packet_source,
-        enable_vpacket_tracking=False,
+        enable_virtual_packet_logging=False,
         enable_rpacket_tracking=False,
         nthreads=1,
         debug_packets=False,
@@ -96,19 +84,9 @@ def __init__(
 
         self.use_gpu = use_gpu
 
-        self.enable_vpacket_tracking = enable_vpacket_tracking
+        self.enable_vpacket_tracking = enable_virtual_packet_logging
         self.enable_rpacket_tracking = enable_rpacket_tracking
 
-        self.virt_packet_last_interaction_type = np.ones(2) * -1
-        self.virt_packet_last_interaction_in_nu = np.ones(2) * -1.0
-        self.virt_packet_last_line_interaction_in_id = np.ones(2) * -1
-        self.virt_packet_last_line_interaction_out_id = np.ones(2) * -1
-        self.virt_packet_last_line_interaction_shell_id = np.ones(2) * -1
-        self.virt_packet_nus = np.ones(2) * -1.0
-        self.virt_packet_energies = np.ones(2) * -1.0
-        self.virt_packet_initial_rs = np.ones(2) * -1.0
-        self.virt_packet_initial_mus = np.ones(2) * -1.0
-
         self.packet_source = packet_source
 
         # Setting up the Tracking array for storing all the RPacketTracker instances
@@ -206,12 +184,13 @@ def run(
             transport_state.geometry_state,
             numba_model,
             transport_state.opacity_state,
-            transport_state.radfield_mc_estimators,
+            transport_state.estimators,
             transport_state.spectrum_frequency.value,
             number_of_vpackets,
             iteration=iteration,
             show_progress_bars=show_progress_bars,
             total_iterations=total_iterations,
+            enable_virtual_packet_logging=self.enable_vpacket_tracking,
         )
 
         transport_state._montecarlo_virtual_luminosity.value[
@@ -268,7 +247,9 @@ def get_line_interaction_id(self, line_interaction_type):
         )
 
     @classmethod
-    def from_config(cls, config, packet_source, enable_vpacket_tracking=False):
+    def from_config(
+        cls, config, packet_source, enable_virtual_packet_logging=False
+    ):
         """
         Create a new MontecarloTransport instance from a Configuration object.
 
@@ -341,9 +322,9 @@ def from_config(cls, config, packet_source, enable_vpacket_tracking=False):
             packet_source=packet_source,
             debug_packets=config.montecarlo.debug_packets,
             logger_buffer=config.montecarlo.logger_buffer,
-            enable_vpacket_tracking=(
+            enable_virtual_packet_logging=(
                 config.spectrum.virtual.virtual_packet_logging
-                | enable_vpacket_tracking
+                | enable_virtual_packet_logging
             ),
             enable_rpacket_tracking=config.montecarlo.tracking.track_rpacket,
             nthreads=config.montecarlo.nthreads,
diff --git a/tardis/montecarlo/montecarlo_numba/base.py b/tardis/montecarlo/montecarlo_numba/base.py
index 12fb8cc8940..987a70fe7b1 100644
--- a/tardis/montecarlo/montecarlo_numba/base.py
+++ b/tardis/montecarlo/montecarlo_numba/base.py
@@ -39,25 +39,33 @@ def montecarlo_main_loop(
     iteration,
     show_progress_bars,
     total_iterations,
+    enable_virtual_packet_logging,
 ):
-    """
-    This is the main loop of the MonteCarlo routine that generates packets
+    """This is the main loop of the MonteCarlo routine that generates packets
     and sends them through the ejecta.
 
+
     Parameters
     ----------
     packet_collection : PacketCollection
-    numba_radial_1d_geometry : NumbaRadial1DGeometry
+        Real packet collection
+    geometry_state : GeometryState
+        Simulation geometry
     numba_model : NumbaModel
     opacity_state : OpacityState
-    estimators : NumbaEstimators
-    spectrum_frequency : astropy.units.Quantity
-        frequency binspas
+    estimators : Estimators
+    spectrum_frequency :  astropy.units.Quantity
+        Frequency bins
     number_of_vpackets : int
         VPackets released per interaction
-    packet_seeds : numpy.array
-    virtual_packet_logging : bool
-        Option to enable virtual packet logging.
+    iteration : int
+        Current iteration number
+    show_progress_bars : bool
+        Display progress bars
+    total_iterations : int
+        Maximum number of iterations
+    enable_virtual_packet_logging : bool
+        Enable virtual packet tracking
     """
     no_of_packets = len(packet_collection.initial_nus)
 
@@ -177,7 +185,7 @@ def montecarlo_main_loop(
     for sub_estimator in estimator_list:
         estimators.increment(sub_estimator)
 
-    if montecarlo_configuration.ENABLE_VPACKET_TRACKING:
+    if enable_virtual_packet_logging:
         vpacket_tracker = consolidate_vpacket_tracker(
             vpacket_collections, spectrum_frequency
         )
diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py
index ba793d79120..d4c21767cbe 100644
--- a/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py
+++ b/tardis/montecarlo/montecarlo_numba/tests/test_cuda_formal_integral.py
@@ -1,27 +1,21 @@
-import pytest
 import numpy as np
-from tardis import constants as c
-from astropy import units as u
-
-from copy import deepcopy
 import numpy.testing as ntest
+import pytest
 from numba import cuda
-from numba import njit
-from tardis.montecarlo.montecarlo_numba import njit_dict_no_parallel
-
 
-import tardis.montecarlo.montecarlo_numba.formal_integral_cuda as formal_integral_cuda
 import tardis.montecarlo.montecarlo_numba.formal_integral as formal_integral_numba
+import tardis.montecarlo.montecarlo_numba.formal_integral_cuda as formal_integral_cuda
+from tardis import constants as c
 from tardis.model.geometry.radial1d import NumbaRadial1DGeometry
-from tardis.montecarlo.montecarlo_numba.numba_interface import NumbaModel
-
-
-from tardis.montecarlo.montecarlo_numba.formal_integral import FormalIntegrator
 from tardis.montecarlo.montecarlo_numba.formal_integral import (
+    FormalIntegrator,
     NumbaFormalIntegrator,
 )
+from tardis.montecarlo.montecarlo_numba.numba_interface import (
+    NumbaModel,
+)
 
-from tardis.montecarlo.base import MontecarloTransportSolver
+from tardis.montecarlo.base import MonteCarloTransportSolver
 
 
 # Test cases must also take into account use of a GPU to run. If there is no GPU then the test cases will fail.
@@ -390,8 +384,8 @@ def test_full_formal_integral(
 
     L_cuda = formal_integrator_cuda.integrator.formal_integral(
         formal_integrator_cuda.simulation_state.t_inner,
-        sim.transport.spectrum.frequency,
-        sim.transport.spectrum.frequency.shape[0],
+        sim.transport.transport_state.spectrum.frequency,
+        sim.transport.transport_state.spectrum.frequency.shape[0],
         att_S_ul_cuda,
         Jred_lu_cuda,
         Jblue_lu_cuda,
@@ -402,8 +396,8 @@ def test_full_formal_integral(
 
     L_numba = formal_integrator_numba.integrator.formal_integral(
         formal_integrator_numba.simulation_state.t_inner,
-        sim.transport.spectrum.frequency,
-        sim.transport.spectrum.frequency.shape[0],
+        sim.transport.transport_state.spectrum.frequency,
+        sim.transport.transport_state.spectrum.frequency.shape[0],
         att_S_ul_numba,
         Jred_lu_numba,
         Jblue_lu_numba,
diff --git a/tardis/montecarlo/montecarlo_transport_state.py b/tardis/montecarlo/montecarlo_transport_state.py
index 4d43f428ee6..8edb5fb7606 100644
--- a/tardis/montecarlo/montecarlo_transport_state.py
+++ b/tardis/montecarlo/montecarlo_transport_state.py
@@ -30,6 +30,18 @@ class MonteCarloTransportState(HDFWriterMixin):
         "last_line_interaction_shell_id",
     ]
 
+    vpacket_hdf_properties = [
+        "virt_packet_nus",
+        "virt_packet_energies",
+        "virt_packet_initial_rs",
+        "virt_packet_initial_mus",
+        "virt_packet_last_interaction_in_nu",
+        "virt_packet_last_interaction_type",
+        "virt_packet_last_line_interaction_in_id",
+        "virt_packet_last_line_interaction_out_id",
+        "virt_packet_last_line_interaction_shell_id",
+    ]
+
     hdf_name = "transport_state"
 
     last_interaction_type = None
@@ -276,12 +288,12 @@ def calculate_reabsorbed_luminosity(
         ].sum()
 
     @property
-    def virtual_packet_nu(self):
-        if self.vpacket_tracker is not None:
+    def virt_packet_nus(self):
+        try:
             return u.Quantity(self.vpacket_tracker.nus, u.Hz)
-        else:
+        except AttributeError:
             warnings.warn(
-                "MontecarloTransport.virtual_packet_nu:"
+                "MontecarloTransport.virt_packet_nus:"
                 "Set 'virtual_packet_logging: True' in the configuration file"
                 "to access this property"
                 "It should be added under 'virtual' property of 'spectrum' property",
@@ -290,12 +302,12 @@ def virtual_packet_nu(self):
             return None
 
     @property
-    def virtual_packet_energy(self):
+    def virt_packet_energies(self):
         try:
-            return u.Quantity(self.virt_packet_energies, u.erg)
+            return u.Quantity(self.vpacket_tracker.energies, u.erg)
         except AttributeError:
             warnings.warn(
-                "MontecarloTransport.virtual_packet_energy:"
+                "MontecarloTransport.virt_packet_energies:"
                 "Set 'virtual_packet_logging: True' in the configuration file"
                 "to access this property"
                 "It should be added under 'virtual' property of 'spectrum' property",
@@ -307,7 +319,7 @@ def virtual_packet_energy(self):
     def virtual_packet_luminosity(self):
         try:
             return (
-                self.virtual_packet_energy
+                self.virt_packet_energies
                 / self.packet_collection.time_of_simulation
             )
         except TypeError:
@@ -319,3 +331,109 @@ def virtual_packet_luminosity(self):
                 UserWarning,
             )
             return None
+
+    @property
+    def virt_packet_initial_rs(self):
+        try:
+            return u.Quantity(self.vpacket_tracker.initial_rs, u.erg)
+        except AttributeError:
+            warnings.warn(
+                "MontecarloTransport.virt_packet_initial_rs:"
+                "Set 'virtual_packet_logging: True' in the configuration file"
+                "to access this property"
+                "It should be added under 'virtual' property of 'spectrum' property",
+                UserWarning,
+            )
+            return None
+
+    @property
+    def virt_packet_initial_mus(self):
+        try:
+            return u.Quantity(self.vpacket_tracker.initial_mus, u.erg)
+        except AttributeError:
+            warnings.warn(
+                "MontecarloTransport.virt_packet_initial_mus:"
+                "Set 'virtual_packet_logging: True' in the configuration file"
+                "to access this property"
+                "It should be added under 'virtual' property of 'spectrum' property",
+                UserWarning,
+            )
+            return None
+
+    @property
+    def virt_packet_last_interaction_in_nu(self):
+        try:
+            return u.Quantity(
+                self.vpacket_tracker.last_interaction_in_nu, u.erg
+            )
+        except AttributeError:
+            warnings.warn(
+                "MontecarloTransport.virt_packet_last_interaction_in_nu:"
+                "Set 'virtual_packet_logging: True' in the configuration file"
+                "to access this property"
+                "It should be added under 'virtual' property of 'spectrum' property",
+                UserWarning,
+            )
+            return None
+
+    @property
+    def virt_packet_last_interaction_type(self):
+        try:
+            return u.Quantity(self.vpacket_tracker.last_interaction_type, u.erg)
+        except AttributeError:
+            warnings.warn(
+                "MontecarloTransport.virt_packet_last_interaction_type:"
+                "Set 'virtual_packet_logging: True' in the configuration file"
+                "to access this property"
+                "It should be added under 'virtual' property of 'spectrum' property",
+                UserWarning,
+            )
+            return None
+
+    @property
+    def virt_packet_last_line_interaction_in_id(self):
+        try:
+            return u.Quantity(
+                self.vpacket_tracker.last_interaction_in_id, u.erg
+            )
+        except AttributeError:
+            warnings.warn(
+                "MontecarloTransport.virt_packet_last_line_interaction_in_id:"
+                "Set 'virtual_packet_logging: True' in the configuration file"
+                "to access this property"
+                "It should be added under 'virtual' property of 'spectrum' property",
+                UserWarning,
+            )
+            return None
+
+    @property
+    def virt_packet_last_line_interaction_out_id(self):
+        try:
+            return u.Quantity(
+                self.vpacket_tracker.last_interaction_out_id, u.erg
+            )
+        except AttributeError:
+            warnings.warn(
+                "MontecarloTransport.virt_packet_last_line_interaction_out_id:"
+                "Set 'virtual_packet_logging: True' in the configuration file"
+                "to access this property"
+                "It should be added under 'virtual' property of 'spectrum' property",
+                UserWarning,
+            )
+            return None
+
+    @property
+    def virt_packet_last_line_interaction_shell_id(self):
+        try:
+            return u.Quantity(
+                self.vpacket_tracker.last_interaction_shell_id, u.erg
+            )
+        except AttributeError:
+            warnings.warn(
+                "MontecarloTransport.virt_packet_last_line_interaction_shell_id:"
+                "Set 'virtual_packet_logging: True' in the configuration file"
+                "to access this property"
+                "It should be added under 'virtual' property of 'spectrum' property",
+                UserWarning,
+            )
+            return None
diff --git a/tardis/montecarlo/packet_source.py b/tardis/montecarlo/packet_source.py
index 5bfd9e64f26..f531a98f59e 100644
--- a/tardis/montecarlo/packet_source.py
+++ b/tardis/montecarlo/packet_source.py
@@ -222,7 +222,12 @@ def create_packet_nus(self, no_of_packets, l_samples=1000):
         xis_prod = np.prod(xis[1:], 0)
         x = ne.evaluate("-log(xis_prod)/l")
 
-        return x * (const.k_B.cgs.value * self.temperature) / const.h.cgs.value
+        if isinstance(self.temperature, u.Quantity):
+            temperature = self.temperature.value
+        else:
+            temperature = self.temperature
+
+        return x * (const.k_B.cgs.value * temperature) / const.h.cgs.value
 
     def create_packet_mus(self, no_of_packets):
         """
diff --git a/tardis/montecarlo/tests/test_base.py b/tardis/montecarlo/tests/test_base.py
index 10ca2c167eb..9c7ec3bffbe 100644
--- a/tardis/montecarlo/tests/test_base.py
+++ b/tardis/montecarlo/tests/test_base.py
@@ -12,30 +12,27 @@
 
 
 @pytest.fixture(scope="module", autouse=True)
-def to_hdf_buffer(hdf_file_path, simulation_verysimple):
-    simulation_verysimple.transport.to_hdf(
+def to_hdf_buffer(hdf_file_path, simulation_verysimple_vpacket_tracking):
+    simulation_verysimple_vpacket_tracking.transport.to_hdf(
         hdf_file_path, name="transport", overwrite=True
     )
-    simulation_verysimple.transport.transport_state.to_hdf(
+    simulation_verysimple_vpacket_tracking.transport.transport_state.to_hdf(
         hdf_file_path, name="transport_state", overwrite=True
     )
 
 
-transport_properties = [
-    "last_interaction_in_nu",
-    "last_interaction_type",
-    "last_line_interaction_in_id",
-    "last_line_interaction_out_id",
-    "last_line_interaction_shell_id",
-]
+transport_properties = [None]
 
 
+@pytest.mark.xfail(reason="No HDF properties being written currently")
 @pytest.mark.parametrize("attr", transport_properties)
-def test_hdf_transport(hdf_file_path, simulation_verysimple, attr):
-    actual = getattr(simulation_verysimple.transport.transport_state, attr)
+def test_hdf_transport(
+    hdf_file_path, simulation_verysimple_vpacket_tracking, attr
+):
+    actual = getattr(simulation_verysimple_vpacket_tracking.transport, attr)
     if hasattr(actual, "cgs"):
         actual = actual.cgs.value
-    path = f"transport_state/{attr}"
+    path = f"transport/{attr}"
     expected = pd.read_hdf(hdf_file_path, path)
     assert_almost_equal(actual, expected.values)
 
@@ -47,12 +44,37 @@ def test_hdf_transport(hdf_file_path, simulation_verysimple, attr):
     "j_estimator",
     "montecarlo_virtual_luminosity",
     "packet_luminosity",
+    # These are nested properties that should be tested differently
+    # "spectrum",
+    # "spectrum_virtual",
+    # "spectrum_reabsorbed",
+    # This is a scalar and should be tested differently
+    # "time_of_simulation",
+    "emitted_packet_mask",
+    "last_interaction_type",
+    "last_interaction_in_nu",
+    "last_line_interaction_out_id",
+    "last_line_interaction_in_id",
+    "last_line_interaction_shell_id",
+    "virt_packet_nus",
+    "virt_packet_energies",
+    "virt_packet_initial_rs",
+    "virt_packet_initial_mus",
+    "virt_packet_last_interaction_in_nu",
+    "virt_packet_last_interaction_type",
+    "virt_packet_last_line_interaction_in_id",
+    "virt_packet_last_line_interaction_out_id",
+    "virt_packet_last_line_interaction_shell_id",
 ]
 
 
 @pytest.mark.parametrize("attr", transport_state_properties)
-def test_hdf_transport_state(hdf_file_path, simulation_verysimple, attr):
-    actual = getattr(simulation_verysimple.transport.transport_state, attr)
+def test_hdf_transport_state(
+    hdf_file_path, simulation_verysimple_vpacket_tracking, attr
+):
+    actual = getattr(
+        simulation_verysimple_vpacket_tracking.transport.transport_state, attr
+    )
     if hasattr(actual, "cgs"):
         actual = actual.cgs.value
     path = f"transport_state/{attr}"
diff --git a/tardis/simulation/base.py b/tardis/simulation/base.py
index 5c817311e1f..40f655824e7 100644
--- a/tardis/simulation/base.py
+++ b/tardis/simulation/base.py
@@ -14,10 +14,11 @@
 from tardis.io.configuration.config_reader import ConfigurationError
 from tardis.io.util import HDFWriterMixin
 from tardis.model import SimulationState
+from tardis.model.parse_input import initialize_packet_source
 from tardis.montecarlo import (
     montecarlo_configuration as montecarlo_configuration,
 )
-from tardis.montecarlo.base import MontecarloTransportSolver
+from tardis.montecarlo.base import MonteCarloTransportSolver
 from tardis.plasma.standard_plasmas import assemble_plasma
 from tardis.util.base import is_notebook
 from tardis.visualization import ConvergencePlots
@@ -369,7 +370,6 @@ def advance_state(self):
         # case it needs some extra kwargs.
 
         estimators = self.transport.transport_state.estimators
-
         if "j_blue_estimator" in self.plasma.outputs_dict:
             update_properties.update(
                 t_inner=next_t_inner,
@@ -695,7 +695,9 @@ def from_config(
                     config, atom_data=atom_data
                 )
             if packet_source is not None:
-                simulation_state.packet_source = packet_source
+                simulation_state.packet_source = initialize_packet_source(
+                    config, simulation_state.geometry, packet_source
+                )
         if "plasma" in kwargs:
             plasma = kwargs["plasma"]
         else:
@@ -711,10 +713,10 @@ def from_config(
                 )
             transport = kwargs["transport"]
         else:
-            transport = MontecarloTransportSolver.from_config(
+            transport = MonteCarloTransportSolver.from_config(
                 config,
                 packet_source=simulation_state.packet_source,
-                enable_vpacket_tracking=virtual_packet_logging,
+                enable_virtual_packet_logging=virtual_packet_logging,
             )
 
         convergence_plots_config_options = [
diff --git a/tardis/tests/test_tardis_full.py b/tardis/tests/test_tardis_full.py
index 55a57370707..5ababcbb6bc 100644
--- a/tardis/tests/test_tardis_full.py
+++ b/tardis/tests/test_tardis_full.py
@@ -97,35 +97,3 @@ def test_virtual_spectrum(self, transport, refdata):
         assert_quantity_allclose(
             transport.transport_state.spectrum_virtual.luminosity, luminosity
         )
-
-    def test_transport_properties(self, transport):
-        """
-        Tests whether a number of transport attributes exist and also verifies
-        their types
-
-        Currently, transport attributes needed to call the model routine to_hdf5
-        are checked.
-        """
-
-        virt_type = np.ndarray
-
-        props_required_by_modeltohdf5 = dict(
-            [
-                ("virt_packet_last_interaction_type", virt_type),
-                ("virt_packet_last_line_interaction_in_id", virt_type),
-                ("virt_packet_last_line_interaction_out_id", virt_type),
-                ("virt_packet_last_line_interaction_shell_id", virt_type),
-                ("virt_packet_last_interaction_in_nu", virt_type),
-                ("virt_packet_nus", virt_type),
-                ("virt_packet_energies", virt_type),
-            ]
-        )
-
-        required_props = props_required_by_modeltohdf5.copy()
-
-        for prop, prop_type in required_props.items():
-            actual = getattr(transport, prop)
-            assert type(actual) == prop_type, (
-                f"wrong type of attribute '{prop}':"
-                f"expected {prop_type}, found {type(actual)}"
-            )
diff --git a/tardis/visualization/tools/sdec_plot.py b/tardis/visualization/tools/sdec_plot.py
index e42814142e8..8490a58b990 100644
--- a/tardis/visualization/tools/sdec_plot.py
+++ b/tardis/visualization/tools/sdec_plot.py
@@ -266,53 +266,55 @@ def from_hdf(cls, hdf_fpath, packets_mode):
             if packets_mode == "virtual":
                 return cls(
                     last_interaction_type=hdf[
-                        "/simulation/transport/virt_packet_last_interaction_type"
+                        "/simulation/transport/transport_state/virt_packet_last_interaction_type"
                     ],
                     last_line_interaction_in_id=hdf[
-                        "/simulation/transport/virt_packet_last_line_interaction_in_id"
+                        "/simulation/transport/transport_state/virt_packet_last_line_interaction_in_id"
                     ],
                     last_line_interaction_out_id=hdf[
-                        "/simulation/transport/virt_packet_last_line_interaction_out_id"
+                        "/simulation/transport/transport_state/virt_packet_last_line_interaction_out_id"
                     ],
                     last_line_interaction_in_nu=u.Quantity(
                         hdf[
-                            "/simulation/transport/virt_packet_last_interaction_in_nu"
+                            "/simulation/transport/transport_state/virt_packet_last_interaction_in_nu"
                         ].to_numpy(),
                         "Hz",
                     ),
                     lines_df=lines_df,
                     packet_nus=u.Quantity(
-                        hdf["/simulation/transport/virt_packet_nus"].to_numpy(),
+                        hdf[
+                            "/simulation/transport/transport_state/virt_packet_nus"
+                        ].to_numpy(),
                         "Hz",
                     ),
                     packet_energies=u.Quantity(
                         hdf[
-                            "/simulation/transport/virt_packet_energies"
+                            "/simulation/transport/transport_state/virt_packet_energies"
                         ].to_numpy(),
                         "erg",
                     ),
                     r_inner=r_inner,
                     spectrum_delta_frequency=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum_virtual/scalars"
+                            "/simulation/transport/transport_state/spectrum_virtual/scalars"
                         ].delta_frequency,
                         "Hz",
                     ),
                     spectrum_frequency_bins=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum_virtual/_frequency"
+                            "/simulation/transport/transport_state/spectrum_virtual/_frequency"
                         ].to_numpy(),
                         "Hz",
                     ),
                     spectrum_luminosity_density_lambda=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum_virtual/luminosity_density_lambda"
+                            "/simulation/transport/transport_state/spectrum_virtual/luminosity_density_lambda"
                         ].to_numpy(),
                         "erg / s cm",  # luminosity_density_lambda is saved in hdf in CGS
                     ).to("erg / s AA"),
                     spectrum_wavelength=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum_virtual/wavelength"
+                            "/simulation/transport/transport_state/spectrum_virtual/wavelength"
                         ].to_numpy(),
                         "cm",  # wavelength is saved in hdf in CGS
                     ).to("AA"),
@@ -322,61 +324,61 @@ def from_hdf(cls, hdf_fpath, packets_mode):
 
             elif packets_mode == "real":
                 emitted_packet_mask = hdf[
-                    "/simulation/transport/emitted_packet_mask"
+                    "/simulation/transport/transport_state/emitted_packet_mask"
                 ].to_numpy()
                 return cls(
                     # First convert series read from hdf to array before masking
                     # to eliminate index info which creates problems otherwise
                     last_interaction_type=hdf[
-                        "/simulation/transport/last_interaction_type"
+                        "/simulation/transport/transport_state/last_interaction_type"
                     ].to_numpy()[emitted_packet_mask],
                     last_line_interaction_in_id=hdf[
-                        "/simulation/transport/last_line_interaction_in_id"
+                        "/simulation/transport/transport_state/last_line_interaction_in_id"
                     ].to_numpy()[emitted_packet_mask],
                     last_line_interaction_out_id=hdf[
-                        "/simulation/transport/last_line_interaction_out_id"
+                        "/simulation/transport/transport_state/last_line_interaction_out_id"
                     ].to_numpy()[emitted_packet_mask],
                     last_line_interaction_in_nu=u.Quantity(
                         hdf[
-                            "/simulation/transport/last_interaction_in_nu"
+                            "/simulation/transport/transport_state/last_interaction_in_nu"
                         ].to_numpy()[emitted_packet_mask],
                         "Hz",
                     ),
                     lines_df=lines_df,
                     packet_nus=u.Quantity(
-                        hdf["/simulation/transport/output_nu"].to_numpy()[
-                            emitted_packet_mask
-                        ],
+                        hdf[
+                            "/simulation/transport/transport_state/output_nu"
+                        ].to_numpy()[emitted_packet_mask],
                         "Hz",
                     ),
                     packet_energies=u.Quantity(
-                        hdf["/simulation/transport/output_energy"].to_numpy()[
-                            emitted_packet_mask
-                        ],
+                        hdf[
+                            "/simulation/transport/transport_state/output_energy"
+                        ].to_numpy()[emitted_packet_mask],
                         "erg",
                     ),
                     r_inner=r_inner,
                     spectrum_delta_frequency=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum/scalars"
+                            "/simulation/transport/transport_state/spectrum/scalars"
                         ].delta_frequency,
                         "Hz",
                     ),
                     spectrum_frequency_bins=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum/_frequency"
+                            "/simulation/transport/transport_state/spectrum/_frequency"
                         ].to_numpy(),
                         "Hz",
                     ),
                     spectrum_luminosity_density_lambda=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum/luminosity_density_lambda"
+                            "/simulation/transport/transport_state/spectrum/luminosity_density_lambda"
                         ].to_numpy(),
                         "erg / s cm",
                     ).to("erg / s AA"),
                     spectrum_wavelength=u.Quantity(
                         hdf[
-                            "/simulation/transport/spectrum/wavelength"
+                            "/simulation/transport/transport_state/spectrum/wavelength"
                         ].to_numpy(),
                         "cm",
                     ).to("AA"),